System and method for estimating a shooter range

[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 missiles, 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 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) from 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. 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. 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 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 invention addresses the disadvantages of known solutions and, in various embodiments, provides methods and systems for estimating shooter range for distant shots, especially when 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.

[0009] The invention is defined in the appended claims.

[0010] According to one aspect of the invention, the method of disambiguating a projectile's trajectory from shockwave-only signals comprises the steps of measuring at least the initial portion of the shockwave-only signals by five or more spaced apart acoustic sensors forming the antenna, estimating the error distribution time of the acoustic sensors, determining from the measured initial portion of the signals only the shockwave of arrival time difference (TDOA) for sensor pairs with a time resolution greater than the estimated time distribution error calculation, and selecting a unified projectile trajectory based on the defined disambiguation confidence interval and the rest value for the TDOA of acoustic sensors.

[0011] According to another aspect of the invention, the method of disambiguating a projectile from shockwave-only signals comprises the steps of measuring at least the initial portion of the shockwave-only signals by five or more spaced apart acoustic sensors forming the antenna, determining from the measured portion of the shockwave-only signals the difference arrival times (TDOA) for sensor pairs, applying a genetic algorithm to the initial chromosome, which involves the assumption of a projectile trajectory, for a fixed number of generations, calculation of residues for solutions obtained with chromosomes from a genetic algorithm, performing a gradient search on a solution having the smallest residue and its ambiguous alternating solution, and if the ratio of the solution having the smallest residue and its ambiguous alternating solution is greater than the specified value, an indication of the solution having the smallest residue as the unified projectile trajectory.

[0012] According to yet another aspect of the invention, in the method of estimating shooter distance by shockwave and muzzle blast detection, shockwave-only signals as well as muzzle blast signals are measured. The initial shooter distance is estimated from the measured shock wave and muzzle blast signals, assuming the initial projectile velocity and the projectile drag coefficient. The instantaneous bullet velocity along the bullet trajectory is calculated iteratively to obtain the updated shooter distance. The number of muzzle blast detection channels is usually smaller than the number of shock wave detection channels.

[0013] Preferred embodiments may include one or more of the following features. Arrival time difference (TDOA) between signals based only on shock wave and muzzle blast signals and the angle of incidence are calculated to determine the initial distance of the shooter. A number of iterations can be made, or the updated shooter distance will be considered the final shooter distance if the relationship between successively determined updated shooter distances meets the convergence criterion. For example, the convergence criterion can be selected in such a way that the difference between the successively determined updated shooter ranges or the percentage change between the successively determined updated shooter ranges is less than a predetermined value. For real solutions, the calculated projectile speed is assumed to always be at least equal to the speed of sound. The solutions are checked for compliance. For example, the updated Sagittarius distance is considered invalid if it turns out that the projectile trajectory angle and the incidence angle are greater than the set value.

[0014] Even if it is determined that the calculated shooter distance is invalid, a solution can still be obtained using a genetic algorithm (GA). For example, you can define an initial population in GA with a fixed number of individuals, in which each individual is represented by a trio that contains the assumed shooter distance, miss azimuth (MA), and miss angle (ME) of the projectile trajectory. GA is performed for a fixed number of generations and the residuals are calculated for individuals in each generation. In each generation, the solution with the smallest residue is chosen as an individual who survives without mutation. The solution with the smallest residue is selected as the updated shooter distance. The solution can be further refined by performing a fixed number of iterations for each trio in a generation to calculate the verified shooter distance, in which the rest for individuals in each generation are calculated using the verified shooter distance.

[0015] GA includes crossover and mutation operators. The crossing operator lists at least one of the missed azimuth and missed angle between two individuals from the population in a given generation, while the mutation operator contains a field mutation (replacing the trio value with a randomly selected value), and an incremental mutation (causing a small mutation in all three fields) , and no mutations (leaving individuals unchanged in a given generation).

[0016] According to yet another aspect of the invention, a method of disambiguating a projectile from shockwave signals and a limited number of muzzle blast signals includes measuring shockwave only signals by five or more spaced-out acoustic sensors, measuring the muzzle blast signals by up to 4 of the sensors, and determining of shockwave-only arrival time difference (TDOA) information for sensor pairs. The method then involves performing a genetic algorithm for a predetermined number of generations with an initial population comprising a specified number of individuals, where each individual is represented by a four that includes shooter azimuth, shooter aim angle, miss azimuth and miss angle, and calculate residuals for individuals in every generation at residuals involving least squares matching of the TDOA combination of shockwave and muzzle blast signals. If the ratio of the solution having the smallest residue and its ambiguous alternating solution is greater than the set value, for example a value of at least 2, the solution having the smallest calculated residue is indicated as the unified projectile trajectory.

[0017] In accordance with another aspect of the invention, a method of extracting a muzzle wave signal in the presence of a shockwave signal comprises defining a time window having a width corresponding to the time necessary for the muzzle wave to pass through the sensor array and detecting the shockwave signal. After detecting the shockwave signal, the window is moved forward in time and the total energy received in the window is measured as a function of the forward time. The maximum measured total energy is associated with the muzzle signal.

[0018] 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 μs . [0019] 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.

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

[0021] In other preferred embodiments, the arrival time difference (TDOA) for the sensor pairs can be determined by designing the sensor that first encounters the shock wave as a reference sensor, and setting the first time trigger 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.

[0022] 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 [0023] 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 have 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 a shockwave basis;

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 Presented [0024] 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 are detected to obtain an accurate solution or when such parameters cannot be reliably detected.

[0025] Trajectories of supersonic missiles are estimated only based on the arrival time of the missile 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 an accurate solution based only on a shock wave is essential for determining the trajectory.

[0026] 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 radius of curvature of the expanding conical surface.

[0027] 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 of 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.

[0028] 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 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 (C_xj, C._sj, C._fromj) on a spherical surface, with one element 22 located in the middle on (Cx0, Cy0, Cz0) so that homogeneous sensor apertures occur irrespective 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. d_and = c ^ (vol_and -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.

[0029] The system preferably has features that ensure that it does 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 using elastomeric shock mounts to isolate it from mast vibration.

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

[0031] Referring now to Fig. 8, the process 800 determines if the detected signal is from a shockwave. Process 800 starts at step 802 and checks at step 804 whether the signal is a loud enough event to be considered an impact, or

the peak signal value exceeds a given parameterized threshold of e.g. 500. If so, process 800 goes to step 806 and checks if there is a sharp transition from zero to the peak signal value, ensuring that the transition to this peak value is not preceded by another signal having a significant magnitude, e.g. 1/16 of the peak signal value.

[0032] 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 if the transition of the peak pressure from 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. around 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.

[0033] Referring again to Fig. 1, it is assumed that the projectile's trajectory coincides with the x axis. 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 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 = n-cos (ę). The angle φ is defined by sin ^) = L / 2 r₁ so that for small angles φ we get φ ~ L / 2 r₁. The time difference between curvature measurement points between the antenna surface, crossing the conical surface with the radius r₁ is equal to dt₁ = Δd / c = (r₁ - d) / c ~ r₁ φ / 2 c = L / (8-yyyy). Difference in curvature measurement time at₂ = year₁ - L ^ sin ^) is given the same ₃ expression, where r₂ replaces r₁. Accordingly, dt₂ = dt₁ + L sin ^) / 8₁²c.

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

[0035] The test measurement made 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 made 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 equal to X is given :

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

(X dt₂)² where p₂(x) =, - e ^2s'

42πσ and p₁(x) (ξ-dh)², 2S² μϊπσ

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

^{P (r}1 < ^r2) = 2

- and

2p - ¥ e ^-at erf (u - a) du dt₁ -DT_{2 =} L³ sin (q) where a =

V2s V28₁²cs [0036] 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 20. Sound speed assumed was c = 340 m / s. 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.

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

[0038] 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, dg V0 where dt = dV / dt [0039] 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.

[0040] 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 angles of error 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 = 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.

[0041] 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 / m.

[0042] 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).

[0043] 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 analog mode, based on the detection of threshold levels, with the digital logic circuits performing 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.

[0044] In another embodiment, the correlation for each channel is calculated 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 [0045] 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).

[0046] The total signal strength is:

-fOO -fOO

Total power = J | h (t) |² dt = JH (f) |² df [0047] 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.

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

[0049] 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:

JI g < ^(T)|^{2 dt}

Total channel power and °

Signal duration [0050] 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_{and, j} (T)

Corr (gi, gj) =

-J (Channel power i) * (Channel power j) [0051] The remainder for the channel i is calculated by calculating the average value for the sensor and after all sensors j:

Rash (i) = average (Σ (1 - Corr (g_t, g_j)) ί &

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.

[0052] 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 10 samples / sec.

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

[0054] In the Cartesian coordinate system with the origin in the center 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}) + V^{M - 1 sin (b}))

Mc where something (b) =

Vx (X0 -Cx) + Vy (X0-Cy) + Vz (X0-Cz)

LMc [0055]

V = ^Vy 'J represents the speed of the supersonic missile

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

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

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

^{10 dt}min = ^{min (vol}obl -^tzm.e ⁾² when the sum is calculated from all sensors.

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

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

[0060] A block diagram of the genetic algorithm (GA) is shown in Figure 7. Instead of starting with one conjecture solution, the GA 70 process 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 one is evaluated 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_mm between the calculated solution and the measured solution.

^dtmin = ^minJΣ ^(tobl -^T: me Ϋ [0061] Intuitively, an algorithm having a low At ™ value is better.

[0062] 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] [0063] [Azimuth, Shot, Aiming Shot] are determined by the angle (θ + β), while [Azimuth, Missing, Angle of Missing] are determined by the angle β (see Fig. Since the muzzle blast is not used in the shock wave-based approach described above, a nominal distance of 100 meters between the sensor array and the shooter is assumed.

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

[0065] 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 taken or further evaluated. Otherwise, step 76 checks if the adopted fitness criteria have been met.

[0066] Fitness criteria may be, e.g., a calculated missed azimuth of <15 ° and / or a ratio of the residues of two synonymous 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.

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

[0068] The following exemplary hybridization and mutation operators were used to demonstrate process 70:

Table 1

<td>Operator Name</td><td>Operator Type</td><td>PROBABLE bieństwo</td><td>Description</td>

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

<td>Crucifixion misses</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>Mutation incremental</td><td>mutations</td><td> 0,4</td><td>Introduces small mutations in all fields chromosome (within &lt;2 ° for shooter information; in within &lt;0.5 ° for misinformation</td>

<td>Mutation switching</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>

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

[0070] In the case of shockwave only detection, the GA algorithm generates a solution using a computer with a 1 GHz clock operating under the Linux operating system 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. miss azimuth &lt;15 ° were ambiguous 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.

[0071] The solutions described above for the projectile trajectory were obtained without muzzle blast detection. 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.

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

[0073] 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 bullet 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 lc0. 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 - 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 its sensor, are written as a function of these geometrical quantities.

A [0074] 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-th sensor, plus. S distance S from the point of propagation to the jth sensor -.

^c0

Al SI D. "I.e. =^t0 +—+— = ^to + - sin^b + 0),

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

[0075] It can be safely assumed that the velocity V of the projectile 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 projectile 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.

[0076] The parameters required for accurate calculation are usually unknown in a real environment, such as e.g. a battlefield. The distance estimation can, however, be significantly improved by using the iterative process 2oo shown in the block diagram in Fig. 1o, which takes into account the decrease in the projectile speed along the trajectory. The 2oo process begins with a 2o2 step under the following assumptions:

co = 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 projectile [0077] 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 speed V0 by Mach number M0, according to the equation

D ^tms ^C0 ⁰ 1 - cos (a) [0078] With 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:

It has _ M0

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

T _ a

V015 [0079] The angle q is associated with a 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) 'tg (b)

JM [0080] 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 _ ^c0 ^(tms + ^Ta) + ^s [0081] 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.

[0082] 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 may be considered valid if distance a traveled by the projectile and distance 5 = a ^{sin (b)} between sin (a) the sensor and the point whose sound wave is radiated from the bullet to the sensor are important numerical values, i.e. not NAN. NAN is a special floating point value representing the result of a computational operation that cannot return a valid numeric value and is usually used to prevent error propagation in the computation process. In addition, both angles α and β should be smaller than the set threshold indicating that the projectile was actually launched towards the sensor matrix.

[0083] As stated above, a couple of numbers (vol_ms, α) is initially used to calculate the shooter distance D0 in zero approximation, bypassing 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 (vol_ms, α), then the solution is rejected.

[0084] 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 measured reliably in most situations. Gunner azimuth and gunner aiming angle, but not shooter distance, can be determined from TDOA for the shockwave only using known sensor matrix coordinates (S_xj, S_sj, S_fromj). It will be assumed that TDOA vol_msbetween the detected shock wave and the muzzle sound can also be measured, so that muzzle sound may not be detected by all sensors.

[0085] 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 and optimization procedures. 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

-t) shock, obl shock, measure + abs

<td>&lt;N</td><td>Γ t</td><td>λ</td><td>λ</td>

<td>Σ</td><td>* ms, calculated</td><td></td><td>-τ</td>

<td>i = 1</td><td>_VM</td><td></td><td>ms, measure</td>

<td>V</td><td>V</td><td>J</td><td>J</td>

[0086] Intuitively, an algorithm that generates smaller At values_mj_n is better.

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

304, where RANGE is the shooter distance D = shown in Fig. 9, MA is the miss azimuth and ME is the miss targeting 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].

[0088] 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 trio with a predetermined azimuth and angle of aim of the shooter, and with the assumed RANGE for a specific trio. 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 MA and ME values ​​for each three, step 308. Distances are also calculated

D = Sj -P between the shooter and each sensor j that detects a shock wave. [0089] For each three, the angle β is calculated by the equation:

cos (b) = A ·^{ρ -} S)

P - S in which the symbol "•" means 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, α and D, step 312. This process iterates is made several times for each of the triples, e.g. three times, as determined in step 312, after which the rest At is calculated for each three_mj_n defined above which includes the muzzle signal, step 314. [0090] 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.

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

[0092] The following crossover and mutation operators were used to demonstrate process 300:

Table 2

<td>Operator Name</td><td>Type operator</td><td>Prawdopod obieństwo</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>Angle crossing aiming misses</td><td>crossing</td><td> 0,5</td><td>Lists the miss targeting angle between two chromosomes</td>

<td>Crucifixion distances misses</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>Mutation incremental</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>

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

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

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

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

[0097] 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 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 the optimization function to disambiguate the solution based only on the shockwave 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.

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

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

dt

Σ (d 'muzzle, calculated on ^{D '} muzzle, measure Γ + Σ ⁽D '

Shock, calculated _{D '}^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).

[0100] 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 with solutions based only on the shock wave, so that the calculated trajectory solution approaches one of the ambiguous solutions more closely , i.e. by making the solutions more clear.

[0101] Without at least some reliable 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.

[0102] For potentially uncertain muzzle blast detection, an attempt can still be made at the beginning to detect muzzle signals, e.g., to 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).

[0103] 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:

<img file="PL2051095T3_D0001.tif" />

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.

[0104] To distinguish it 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.

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

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