Systems and methods for disambiguating shooter locations

[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 shooters 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 missiles, such as small arms and artillery shells, by measuring parameters associated with the shock wave produced by the projectile. One such system, as described in US Patent No. 5,241,518, consists of at least three spaced sensors, each sensor having three acoustic transducers arranged in a plane. Sensors generate signals in response to the shock wave, which are related to the azimuth and angle of aiming and to the start point of the shock wave. 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 aim of the shooter with respect to the position of the sensor is usually determined by measuring the time of arrival (TOA) information 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) of the sensor (sensors) in relation 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., 1 m away), or very dispersed (e.g. mounted on a vehicle, or carried by soldiers on the battlefield), and can measure outlet pressure and shockwave in an omnidirectional manner in their respective positions. Nevertheless, if the sensors are not relatively wide 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 transformed to map 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 greater than about 300 m.

[0007] US 5 930 202 discloses a method for determining the unambiguous trajectory of a projectile through a matrix of spaced-up acoustic sensors according to the preamble of claim 1.

Summary of the Invention [0008] The invention addresses the disadvantages of known solutions, providing a method characterized according to claim 1. 1.

[0009] Preferably, the method further comprises performing a predetermined number of generations of the genetic algorithm with an initial population that includes the predetermined number of individuals, each individual representing a four that includes shooter azimuth, shooter aiming angle, miss azimuth and miss targeting angle, and calculation of residuals for individuals in each generation, these residues contain least squares matching of the combination of shock wave TDOA measurements and muzzle blast signals. If the ratio of the solution having the smallest residue and its ambiguous alternative solution is greater than the set value, for example a value of at least 2, then the solution having the smallest calculated residue as the unambiguous projectile trajectory is determined. 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 has positive values ​​is less than the minimum time or greater than the maximum time, e.g. less than about 70 μs or greater than about 300 μs.

[0010] 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 maximum total energy, and if the peak signal value is greater than the measured total energy range by a fixed ratio, the peak signal value can be identified as associated with the muzzle signal.

[0011] 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.

[0012] In other preferred embodiments, the time difference of arrival (TDOA) for sensor pairs can be determined by designing the sensor that first encounters the shock wave as a reference sensor, and setting the first time-measuring system flip-flop when the amplitude of e.g. the initial part of the signal only 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 in relation to the reference sensor are then recorded.

[0013] 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 [0014] 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;

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

Fig. 3 schematically shows the ambiguity inherent in determining trajectories based on shockwaves only; Figure 4 schematically shows the probability density for measuring arrival time difference to determine Mach cone curvature; Fig. 5 schematically shows the probability of correct unification between shooter trajectories;

Figure 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;

Figure 8 is a block diagram of the process of discriminating between signals not related to a shockwave;

Figure 9 schematically shows a shock wave arrival time (TOA) model; Figure 10 is a block diagram of a distance estimation process; Figure 11 is a block diagram of a genetic algorithm process for distance estimation.

Detailed Description of Some Embodiments Presented [0015] As described above in the Summary, the invention relates, in various forms, 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 are detected to obtain an accurate solution, or when such parameters cannot be reliably detected.

[0016] 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.

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

[0018] The angle of incidence of the conical surface that first reaches the antenna defines two possible relative angles (often referred to as "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 shooter angular positions.

[0019] The required accuracy can be estimated by considering the shock wave propagation characteristics of Figure 1. Referring now to Figure 2, the antenna 20 includes N sensors (N = 7) capable of determining the arrival times of a progressive cone shock wave. Since it can basically be expected that the trajectories of an incoming projectile can take their origin anywhere, the 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 time 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 successive moments marked as i.e. Distance of sound propagation in the direction of the progressing conical surface is obtained by multiplying each of the time differences by the local sound speed c, i.e. dj = c ^ (t_j -this). 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 stated above, errors in arrival time and sensor coordinate measurements 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.

[0020] 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 electronic circuits, which can also be attached to the sensor mast with elastomeric shock mounts to isolate it from mast vibration.

[0021] 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.

[0022] Referring now to Fig. 8, the process 800 determines whether the detected signal is from a shockwave. The process 800 starts from 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 for 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 value.

[0023] If this is the case, the process 800 proceeds to step 808 and checks if the time between the shockwave minima and maxima is of a sufficiently large value, 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 that 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.

Referring again to Figure 1, it is assumed that the projectile's trajectory coincides with the x axis. The Mach angle is given by the equation dany = arcsin (1 / M), where M is the Mach number determined as the projectile speed V divided by the speed of sound c. L is the characteristic length of the antenna. The radii of curvature of the cone at the two ends of the antenna 20 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 = r ^ cos (f). The angle f is defined by sin (f) = L / 2 r1 so that for small angles f we get f ~ L / 2 r1. 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₁ f / 2 c = L / (8 »ri» c).

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

[0025] 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. Exemplary distributions of measurement values ​​for the time differences dt1 and dt2 are shown in Fig. 4.

[0026] 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}and < ^r2 ^{1 x)} = P2 (^x) JPL ^{(X) dX} x

where p₂(x) =

42Ps (xd / e ^2S and p₁(x) (X-dti)², 2σ²

V2Ps [0027] Integration after x and substitution of variables results in the following probability of making the correct decision:

1 ^¥_-at 2

P (r <r₂) = --- j = I e ^at erf (u -a) big π -_¥ dt₁ - dt ₂ L³ sin (q) where a = ' ² = ₂V2s y / 28₁cs [0028] 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 the antenna A sound speed 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 high CPA values, the curvature difference at the two antenna ends (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.

[0029] 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.

[0030] 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, at dg V0 dt dV / dt [0031] 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 replacing the relative change in gain 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.

[0032] 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 ^-3

340 ·> / Ϊ2 = 0.85ms. Standard deviation of 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 = 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.

[0033] 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 ^ / θ, 85² + 0,75² + 0,25² = 1.1 / m.

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

[0035] In one embodiment, the exemplary system uses an analog time difference of arrival (TDOA) system that uses 1 MHz clocks on 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 analogue 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 the next shockwave.

[0036] In another embodiment, a correlation is calculated for each channel with each other channel for a time period centered on the time of TDOA hardware detection. 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 [0037] 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).

[0038] The total signal strength is:

+¥ +¥

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

[0040] 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.

[0041] Referring now to Fig. 6, in process 60, shockwave time signal data gi (t), gj (t) is 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 [0042] 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_{i, j}(T)

C_orr (g, g) =. ^j^{1 1} yJ (Channel power i) * (Channel power j) [0043] The remainder for channel i is calculated by calculating the average value for the sensor and after all sensors j:

Re (i) = average £ (1 - Corr (g_t, g_})) as specified in step 612. TDOA and correlations for each channel with the best (ie least 12) integer residue is then selected as the "best" solution, step 614.

[0044] As said above, channel data is typically sampled at discrete time intervals with a predetermined sampling frequency of, e.g., 41.666.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 / sec.

[0045] Once the individual arrival time differences (TDOA) 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. [0046] In the Cartesian coordinate system with the origin in the center of the matrix, i.e. {(Cx0, Cy0, Cz0) = (0, 0, 0)}, the time of arrival of the shock wave TOA for a given sensor (Cxj, Cyj, Cyj) (see Fig. 2) is given by:

^t Shock ^{= t}0 ⁺

L

Mc (cos (b) + VM² -1 sin (b)) with cos (b) (V Ί

Vx (X0 - C) + Vy (X0 - Cy) + Vz (X0 - C) LMc ^Vs

J represents the speed of the supersonic missile

Mc = V = yjVj + Vy + Vj, where c is the speed of sound and M is the Mach number, β represents the "box angle" between the shooter's position and the projectile's trajectory, which includes both azimuth and aiming angles. A direct hit will be β = 0. The Mach angle q defines the relation 1 / M = sin (Q).

[0047] 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 TDOA resolution is high enough to solve the shockwave curvature, then these two almost identical solutions can be disambiguated. The fundamental ambiguity of TDOA solutions based only on shockwaves is shown in Figure 3.

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

Dtn = ^min ^ Σ ^(toBL ^-tmeasure when the sum is calculated from all sensors.

[0049] One approach to solve this problem is the Levenberg Marquardt L1 algorithm, described in detail in US Patent No. 5930202. 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 search direction based on a pre-determined transition rule, such as direct methods using the objective function and limit values, and gradient methods using derivatives of the first and second degree. 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.

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

[0051] 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 the Euclidean distance Atmi_n between the calculated solution and the measured solution.

^dt = ^min ^ Σ ^(toBL ^-tmeasure [0052] Intuitively, an algorithm having a small value of At ,,,,,, is better.

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

[Azimuth Shot, Aiming Angle, Shot, Azimuth, Missing, Aiming Angle] [0054] [Azimuth, Shot, Aiming Shot] are determined 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 meters between the sensor matrix and the shooter is assumed.

[0055] 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}.

[0056] In step 75, it is checked that the achievement is the maximum number of iterations for GA, set eg to 25. 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 estimated. Otherwise, step 76 checks if the adopted fitness criteria have been met.

[0057] Fitness criteria may be, e.g., a calculated missed azimuth of <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.

[0058] In each generation, the "best" individual can survive unchanged, while the 100 best individuals, assessed according to their fitness criteria, also survive, but they are used to create another 100 individuals from pairs of survivors with the help of crossover operators / mutations listed in Table 1.

[0059] 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 the shooter trajectory / azimuth between two chromosomes</td>

<td>Missing crosses</td><td>crossing</td><td> 0,5</td><td>Lists 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>

[0060] Disambiguation is obtained and / or improved by performing gradient search in the best solution and in a suitable alternative solution. Residues and residue ratios are calculated for both ambiguous solutions. If the calculated azimuth of miss 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. [0061] In the case of shockwave detection only, the GA algorithm generates a solution using a computer with a 1 GHz clock running 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. those with 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 geometry of the sensor matrix and the assumed distribution of shots, with shots having a small aiming angle easier to disambiguate.

[0062] The solutions described above for the projectile trajectory were obtained without detecting the muzzle blow. 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.

[0063] Fig. 9 schematically shows the arrival time model (TOA), which is described in more detail in US Patent No. 6,178141. The TOA model can be used to estimate the projectile trajectory and shooter direction relative to the sensor location. 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.

[0064] 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), individual sensors are located at S points _j (Sxj, Syj, Szj), and the projectile trajectory is shown as running from the shooter in the direction of A. The vector from the shooter to the jth sensor is marked with 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 = V / c0. M is the mach number of the projectile, V is the supersonic speed of the projectile, and c0 is the speed of sound (depending 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 20 y plane, and the aiming angle measured up from the x - y plane. The equations defining the time of arrival of the shock wave, i.e. the unit vector at the j-th sensor, were written as a function of these geometrical quantities.

iai [0065] The arrival time is equal to the time it takes the bullet to travel the distance A to the point from which the sound propagated towards the jth sensor, plus the distance S from the propagation point to the jth sensor.

Al SI Dl. "I.e. =^t0 + —+ —= ^t0 + ^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.

[0066] It can be safely assumed that the velocity V of the bullet remains constant on the path corresponding to the sensor spacing, such that there is a slight decrease in speed after 17 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.

[0067] 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

M ₀ = V₀ / c = initial Mach number of the bullet [0068] The first estimation of the shooter distance Do 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 array and assuming an initial constant velocity V0 by Mach number M0, according to the equation fns' ^C0

D0 =

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

c ^ lv_0} so that the time the bullet travels the distance along the trajectory, step 208, is expressed as:

T =

V0 [0070] The angle q is related to the Mach number Ma with the equation:

^{sin (Q)} = 77” ^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) · tg (b)

X² [0071] 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 Ma and updated flight time Ta, for distance a. The measured TDOA tms and the calculated updated values ​​Ta and a are then used for the next update of the value D of shooter distance:

D = c0 (tms + Ta) + s [0072] 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.

[0073] Process 200 also checks in step 214 whether the updated D value of the distance between the shooter and the sensor matrix is ​​a "reasonable" value, in which case the process 200 ends in step 216. For example, the value of D can be considered valid if the distance traveled by the projectile and distance 5 = 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 representing the result of a calculation operation that cannot return a valid numerical 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.

[0074] As stated above, a pair of numbers (tms, α) are initially used to calculate the shooter distance D0 in zero approximation, ignoring the change in projectile speed 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.

[0075] 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.

[0076] 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 -t^j )² shock, obl shock, measure —-- + abs (n (t 7 ms, calculated

-t ms, measure [0077] Intuitively, the algorithm generating smaller values ​​of Δτ ^ is better.

[0078] 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 exemplary 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].

[0079] The calculations are carried out according to a process similar to the one outlined earlier with respect to the solution based only on the shock wave. 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.

[0080] 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 over 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, a, 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 remainder Atmin defined above which contains the muzzle signal is calculated for each three, step 314.

[0081] 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.

[0082] In each generation, the "best" individual can survive unchanged, while the 100 best individuals, assessed according to their fitness criteria, also survive, but they are used to create another 100 individuals from pairs of survivors with the help of crossover operators / mutations listed in Table 2 below. [0083] 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 aiming angle miss 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 field) randomly selected new value from 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>

[0084] The GA 300 process is performed with an initial population of 200 different triples, with a degree of completion 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.

[0085] 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 shock wave.

[0086] Calculations that do not take into account the decrease in projectile speed along its path 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.

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

[0088] 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 a modified cost function or the rest, so that any available muzzle signals are "mixed in" to the optimization function to disambiguate the solution based only on the shock wave and / or to refine the estimation of the shooter distance. Nevertheless, if less than 3 muzzle and less than 5 blast channels are detected, an alarm can be triggered without attempting to locate the shooter.

[0089] 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 means that the system binds it to "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.

[0090] If only muzzle blast signal is detected on some channels, the rest in this case can be defined as ^dtmn = Σ ^(Dtmuzzle, calculated ^{- Dt}mzzZe, measure,) + Σ Σ ^ \^dtShock calculated ^dt 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).

[0091] 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 solution of the trajectory approaches one of the ambiguous ones more closely solutions, i.e. making the solutions more unambiguous.

[0092] 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.

[0093] 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 (offset by a certain number of frequency intervals to account for delays filters).

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

· Look 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:

<img file="PL2199817T3_D0001.tif" />

n = 0 j = 0 where square mb_andⁿ+ j represents the amount of energy, e.g. muzzle blast energy, measured as 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.

[0095] 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 certain 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.

[0096] The detected muzzle blast signal can then be used to determine the shooter distance and / or disambiguate the shockwave signal as described above. [0097] 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.

[0098] 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. [0099] Although the invention has been disclosed in connection with preferred embodiments shown and described in detail, various changes and improvements can be made thereto. As an example, although the characters depicted are shown to have acoustic sensors such as microphones, this need not be the case. Instead, other types of mechanical or electrical pressure sensors can be used. Moreover, the values ​​given in Tables 1 and 2 for different operators are intended only as an example, while other values ​​can be chosen depending on the actual field conditions.