Minimum miss distance vector measuring system
This is a Continuation-in-Part of application Ser. No. 565,514 filed on Apr. 7, 1975 and now abandoned. 
    
    
     FIELD OF THE INVENTION 
     The invention relates to the solution of the minimum miss distance vector problem of a missile trajectory past a target. 
     BACKGROUND OF THE INVENTION 
     Matrix inversion techniques utilized in attempts to solve this problem before have yielded poor results because of the instability of the mathematical model and the indeterminate nature of the required matrix inversion. Simple triangulation methods fail because of discontinuities and insufficient baseline lengths to provide adequate accuracy. A Gauss-Newton estimation procedure has been the classic approach to the problem. (See Ortega and Rheinboldt, infra, at p. 267.) This incorporates the use of quasi-linear estimation techniques. The short baseline lengths of these systems account for excessive sensitivity of trajectory parameters on very small range errors, manifested in a highly ill-conditioned covariance matrix, making the estimate inaccessible. 
     SUMMARY OF THE INVENTION 
     In the present invention, space diverse electronic range sensors are mounted on the target to sequentially sense a plurality of ranges to the missile as the missile approaches the target. The range data so accumulated is communicated, together with corresponding sensor identification data, to a digital data processor. The data processor, utilizing one of several possible nonlinear estimation algorithms, iteratively establishes the minimum miss distance vector of the missile with great accuracy and in a relatively short period of time. Any of the mathematical methods allows for missing data and is very stable in operation. 
     It is an object of the invention to provide a plurality of range and time data points for a missile having a trajectory in the vicinity of a target. 
     It is a further object of the invention to utilize digital data processing techniques to derive a minimum miss distance vector for the missile. 
     It is still a further object of the invention to provide the minimum miss distance vector in a short period of time, in the absence of some data pairs and in a stable manner. 
     It is an additional object of the invention to provide the minimum miss distance vector in the absence of continuous or unambiguous data. 
    
    
     DESCRIPTION OF THE DRAWINGS 
     FIG. 1 illustrates the operational configuration of the system of the invention. 
     FIG. 2 illustrates the operational block diagram of the &#34;steepest descent&#34; algorithm which may be utilized in the data processor of the invention. 
     FIG. 3 illustrates the operational block diagram of the &#34;N-step Conjugate Gradients&#34; algorithm which may be utilized in the data processor of the invention. 
     FIG. 4 illustrates the essential detection geometry of the system of the invention. 
     FIG. 5 illustrates in a more detailed block diagram form, synchronizer 28 of FIG. 1. 
     Telemetry Data Handling Unit 68 of FIG. 5 may be designed by one having average skill in the art. 
    
    
     DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT 
     Referring to the drawing, it will be seen that the target aircraft 10 has four antennas 12, 14, 16 and 18 mounted respectively on a rudder tip, each of the wing tips and a point well forward. These antennas are fed by Receiver-Transmitters (R/T units) 20, 22, 24 and 26, respectively, of a multi-signal radar system. Each of the R/T units emits a radar pulse in sequence. These pulses may be, for example, 40 nanoseconds long and they are transmitted sequentially under control of synchronizer 28 on the target which may be an aircraft. The time spacing between the radar pulses may be, for example, 400 nanoseconds. This allows reflections from missile 30 to return to the receiver associated with the emitting transmitter before the next transmitter pulse is emitted. This is true because the maximum distance range necessary may be on the order of 185 feet. 
     The synchronizer 28 may incorporate range gate functions to prohibit reception of range signals other than from desired ranges. In the preferred embodiment of the invention, the maximum range is limited to 185 feet and the range gates are programmed to step in increments of 1/4 foot. 
     Time separation of the four R/T unit pulses avoids the necessity for operating the units on different frequencies or otherwise identifying a particular return pulse with a given transmitted pulse. The utilization of the range gate stepping system also avoids excessive extraneous noise in the system. 
     As a pulse is emitted from each one of the R/T units 20, 22, 24 and 26, the synchronizer 28 starts a range counter. The pulse signal is reflected 34 from the missile and received in the same R/T unit where it is converted to a video signal. Each of these video signals is then fed to the synchronizer 28 and each is used to stop a range counter, thereby creating a digital signal which is proportional to the range between the missile 30 and the target 10. 
     An appropriate digital word is generated in synchronizer 28, incorporating this digital range and a digital code which serves to identify the particular R/T unit and antenna from which the range data was derived. 
     Referring to FIG. 5, PRF oscillator 60 generates a frequency of, for example, 1.6 megahertz. This frequency is divided by four in circuit 62. Each of sensors 20, 22, 24 and 26 (FIG. 1) is synchronized to transmit upon receipt of each fourth oscillator pulse from PRF oscillator 60. Divide by four circuit 62 furnishes a two bit code which (a) identifies which of sensors 20, 22, 24 or 26 triggers range counter 64, (b) allows steering of range gate pulses by sensor steering unit 66 to the appropriate sensor, and (c) allows telemetry data handling unit 68 to tag each range detection with the appropriate sensor identification. 
     FIG. 5 synchronizer 28, also shows sensor 20 (See FIG. 1). It will be understood that sensors 22, 24 and 26 operate in the same manner as sensor 20. Divide by four circuit 62 outputs a two bit PRF code to transmitter 70. It will be understood that this two bit code may comprise four different combinations, 00, 01, 10, and 11, on successive PRF input pulses. Transmitter 70, for example, may respond to, for instance, the 00 code. Transmitters in the other three sensors 22, 24 and 26 will, of course, each respond to one of the other three two bit code combinations. When transmitter 70 recognizes, for example, the 00 code, it is enabled to output a transmitter pulse at the frequency of frequency source 72 through circulator 74 to antenna 12. Transmitter 70 accomplishes this output by gating a portion (approximately 40 nanoseconds) of continuously running frequency source 72 to antenna 12. This radio frequency pulse is transmitted 32 to missile 30 (see FIG. 1) and returned 34 back to antenna 12. Antenna 12 feeds this signal through circulator 74 to receiver 76. Receiver 76 amplifies the signal and mixes it with a sample of transmitted frequency from frequency source 72 by means of coupler 78. Circulator 74 serves as a duplexer connecting transmitter 70, antenna 12 and receiver 76 in the proper relationship, as is well known in the art. The output of receiver 76 is a series of bipolar video pulses at the PRF rate and with a width commensurate with the transmitter gate width, for approximately 40 nanoseconds as above-mentioned. These pulses are fed to N range gate channels 80, 80&#39; in the signal processing section of the sensor. It will be understood that the number of N-range gate channels will be determined by the required accuracy and maximum range of the system. 
     A portion of the transmitter pulse will be amplified by receiver 76, converted to a video pulse and be fed through sensor selector 84 to trigger monostable 86. The duration of monostable 86 is slightly longer than the time required for a reflected signal to be received from maximum range. Monostable 86 then opens gate 88 allowing counter 64 to count cycles of range oscillator 92, which may be at a frequency of, for example, 250 megahertz. A number of sample times (N) are generated by decoding counter output 94. At each desired time (corresponding to a range from sensor 20), monostables 96, 96&#39; are triggered, forming a sampling pulse approximately 40 nanoseconds wide. This pulse allows a range gate to sample receiver 76 output at a fixed range for many PRF intervals, thus recovering pulse amplitude modulation at a Doppler frequency rate if target 30 (see FIG. 1) is present at the selected range. Range counter 64 is reset to zero each PRF interval by output 98 from PRF oscillator 60. 
     Target detection for each range interval is accomplished by feeding output 79 of receiver 76 through a signal processing chain comprising range gates 80, 80&#39;, Doppler filters 100, 100&#39;, detectors 102, 102&#39;, low pass filters 104, 104&#39; and threshold devices 106, 106&#39;. The number, N, of signal processing channels is determined by the accuracy and maximum range desired. Range gates 80, 80&#39; (sometimes referred to as boxcar circuits) recover an audio signal generated by Doppler shift of moving target 30 (see FIG. 1). This signal is filtered, detected (rectified) and fed to low pass filters 104, 104&#39; with a time constant much lower than the period of the lowest Doppler frequency expected. Presence of the target in the prescribed range interval will ultimately allow output of low pass filters 104, 104&#39; to exceed a threshold, allowing indication of target presence to telemetry data handling unit 68. In addition to multiplexing data for serial transmission through a telemetry link to the ground station telemetry data handling unit 68 may assign a time tag to each range detection. 
     Since the R/T units 20, 22, 24 and 26 sequentially sense the ranges to the missile, a digital counter 64 provides data to synchronizer 28 corresponding to each R/T unit. The data from counter 64, identified as to which R/T unit and antenna it was derived, is then used to digitally modulate a carrier signal used to transmit the data pairs to a remotely located data processor 36. (See FIG. 1.) 
     Of course, data processor 36 does not have to be located remotely, but could be located in, on, or near target 10. In these cases, wire connections may be used to connect the data output from synchronizer 28 to data processor 36. 
     However, in the case of remote location of the data processor 36, the digital signal is demodulated at the remote location and fed to data processor 36. Data processor 36 adds time of acquisition data to each segment of range-sensor identification data received. The data processor is also provided with information as to the relative positions of the antennas on the target vehicles. 
     Data processor 36, part of ground station 38 is programmed to provide a mathematical solution for trajectory 40 of missile 30 with respect to target 10 and to provide the minimum miss distance vector of missile trajectory 40. 
     One of several mathematical methods known as &#34;Conjugate Directions&#34; may be utilized to accomplish, by an iterative process, the solution of the miss distance vector problem in conjunction with the system of the invention described herein. 
     The first method is the one commonly known in the art as the &#34;Steepest Descent&#34; Algorithm. This algorithm is well known in the art and, for example, may be found completely described by Ortega, J. M. and Rheinboldt W. C., Iterative Solution of Nonlinear Equations in Several Variables, Academic Press, 1970, p. 245. The second method is one commonly known as the Conjugate Gradients Algorithm. This algorithm is also well known in the art and, for example, may be found completely described by Hestenes, M. R. and Stiefel, E., &#34;Methods of Conjugate Gradients for Solving Linear Systems&#34;, Journal of Research of the National Bureau of Standards, 1952, Vol. 49, No. 6, pp. 409-436. Either method involves estimations of the trajectory by dynamic triangulation means followed by direct computation of the minimum miss distance vector. 
     For all practical purposes it suffices to assume a quadratic path model for the relative missile trajectory, namely, 
     
         p(t) = at.sup.2 + vt + s                                   (1) 
    
     where p(t) represents the 3-dimensional relative trajectory vector, while a, v and s stand for three-dimensional relative acceleration, velocity and displacement vectors comprising the nine-dimensional trajectory state vector ##EQU1## The paragraphs which follow are concerned with the fundamental problem of estimating x from measured range and time data and the subsequent determination of the minimum miss distance vector. 
     It is shown below that the estimation of x is formulated as a minimum-seeking problem. The ensemble of measured data is incorporated into a functional F(x) having a minimum value at an optimal estimate x corresponding to the least-squares solution. 
     Assume a total of N detections. With reference to FIG. 2 the ith detection simply states that 
     
         ∥R.sup.i ∥.sup.2 - R.sub.i.sup.2 = 0 ; i = 1 . . . N (3) 
    
     where 
     r i  = p(x, T i ) - α i   
     r i  = the range vector at t = T i   
     p(x, T i ) = trajectory vector at t = T i   
     α i  = antenna position vector at ith detection 
     R i  = scalar range at t = T i   
     Unlike a linear system of equations, the quadratic system (3) does not lend itself to direct root-finding methods. Instead, one may get a least-squares approximation by simply solving an &#34;equivalent&#34; minimum seeking problem involving the minimization of a functional associated with system (3). 
     A convenient functional is derived as follows. Corresponding to the set of N detections, define the error functions 
     
         F.sub.i (x) = ∥r.sup.i ∥ .sup.2 -R.sub.1.sup.2 ; i = 1 . . . N                                                     (4) 
    
     construct a functional by forming some convenient combination of these functions whose minimum constitutes a compromise to minimizing each F i  individually. One such functional is: ##EQU2## namely, the unweighted sum of squared error functions. The value of x that minimizes F (x) constitutes a least-squares approximation to system (3). A weighted functional of the form ##EQU3## where: W i  = weighting coefficient 
     may be used, allowing for stochastic nonlinear optimal estimation. 
     The weighting coefficient, W i , is given by: 
     
         W.sub.i = 1/E {[∥r.sup.i ∥ .sup.2 - (R.sub.i + ξ.sub.i).sup.2 ].sup.2 } 
    
     where: 
     ξ i  = random variable representing ith measurement error. 
     R i  = exact scalar range at t = t i . 
     E = the expectation operator 
     Note that R i  + ξ i  = R i , the measured scalar range at T = T i . When ξ i  is assumed to be a normally distributed random variable having a mean μ i  and variance σ i   2 , the ith weighting coefficient may be shown to be, specifically, 
     
         W.sub.i = 1/[4R.sub.i.sup.2 (σ.sub.i.sup.2 +μ.sub.i.sup.2) + 4R.sub.i μ.sub.i (3σ.sub.i.sup.2 +μ.sub.i.sup.2) + 3σ.sup.4 +6σ.sub.i.sup.2 μ.sub.i.sup.2 +μ.sub.i.sup.4 ] 
    
     The analysis which follows applies to equation (5), above. If the analysis is to be applied to equation (6), above, W i , the weighting coefficient must be added as a multiplying factor within the summation of each of equations (7), (8), (12), (13), (14) and (15), below. 
     Following is a description of three numerically-stable parameter optimization procedures useful for minimizing F(x) and generating an optimal estimate x that characterizes the missile trajectory relative to the target. 
     The three optimization methods discussed below are classified as descent or relaxation methods which start with an initial guess for x and subsequently generate improved estimates by optimally relaxing F(x) along intrinsic search directions in an iterative manner, eventually producing an estimate sufficiently indistinguishable from the optimal solution. 
     The Steepest Descent method is the simplest of the three parameter optimization techniques considered. It is characterized by optimal relaxation along negative gradient directions. [The gradient of a scalar function F(x) is the vector of partials ∇ x  F(x) pointing in the direction of maximum increase of F(x) from point x. As such, the gradient represents the sensitivity of F(x) with respect to x.] 
     The specific algorithm for minimizing F(x) is given below: 
     i. Given estimate x 
     ii. Compute gradient vector ##EQU4## iii. Compute optimal step-size λ from ##EQU5## iv. Update current x by 
     
         x = x - λ γ 
    
      and return to (ii). 
     The algorithm is repeated until γ has reached a sufficiently small neighborhood of zero whence subsequent iterations do not add discernably to x. 
     A visual aid toward understanding the filtering process of the algorithm is given in FIG. 2 in the form of its functional block diagram. Input m represents the measurement vector; in this case, the Sensor-Range-Time data. Input x stands for the current estimate of the state vector. The gradient generator simply takes m and x and produces the gradient or sensitivity vector γ. A two-way switch first presents γ into a step-size generator, which along with x produces the optimal step-size λ [may be thought of as the optimal gain of the feedback amplifier] which, in turn, multiplies the subsequently switched γ resulting in the updating step λ γ. The current estimate x is now updated to x - λγ by means of the update loop in a manner regulated by the three-way switch there. Included in the block diagram are two convergence indicators, namely, the functional value F(x) and the gradient magnitude ∥ γ ∥. Note that, unlike a common feedback controller, the Steepest Descent controller employs a feedforward loop that presents x into the step-size generator; without it, λ could not be determined nor could stability be guaranteed. 
     The reader is invited to turn his attention to the actual computations needed to implement the Steepest Descent process. It can be shown that for the functional: ##EQU6## the gradient vector is given by ##EQU7## in view of (4) ##EQU8## where ##EQU9## Combining in (8) yields ##EQU10## In compact Kronecker notation (12) takes the alternate form ##EQU11## where ##EQU12## and   is the symbol denoting Kronecker multiplication of vectors. (See, Bellman, R., Introduction to Matrix Analysis, McGraw Hill, 1960, pp. 223-239.) The optimal step-size is simply that value of λ which minimizes F(x - λγ). This is a single-variable minimization problem carried out as indicated below. Explicitly, ##EQU13## where: A i  = ∥M i  γ∥ 2   
     B i  = 2r i .spsp.T 
     c i  = F i . 
     evidently a quartic function of the parameter λ. To get the optimal value of λ, simply set the first λ-derivative to zero and solve the resulting equation, namely ##STR1## Using the Newton&#39;s method (see Ortega and Rheinboldt, supra) this can be solved for the three possible roots numerically. The root appropriate for this purpose is the smallest real positive root. This choice insures optimal descent within a convex neighborhood of the search. 
     Although the Steepest Descent method is numerically stable, it is by no means efficient in the sense of convergence speed. In contrast with the Steepest Descent method, other more sophisticated parameter optimization techniques are known to guarantee convergence within a finite number of iterations. In their paper, M. R. Hestenes and E. Steifel supra, introduce the method of Conjugate Gradients and show that convergence may be attained within a number of iterations not exceeding the dimensionality of x, provided that F(x) is a quadratic functional. The corollary here is that only a finite number of iterations are needed for the quartic F(x) in the present case. 
     The specific Conjugate Gradients algorithm for minimizing a given functional F(x) with respect to x is as follows: 
     i. Given estimate x 
     ii. Compute gradient vector ##EQU14## iii. Using the previous gradient vector γ, update the current search vector by ##EQU15## iv. Compute optimal step-size λ from ##EQU16## v. Update current estimate by 
     
         x = x + λs 
    
      and return to (ii). 
     A functional block diagram for the Conjugate Gradients process is given in FIG. 3. 
     One variation of the Conjugate Gradients algorithm involves a somewhat simpler updating formula for the search vector s, namely, ##EQU17## resulting in the so-called One-step Conjugate Gradients method. 
     As its name might imply, the N-Step Switched Conjugate Gradients Method variation of the Conjugate Gradients Method consists of using the normal update formula for s throughout blocks of N consecutive iterations, at the end of which s is reset to zero. This scheme is numerically efficient. 
     With the relative trajectory vector estimate x at hand, it is now possible to determine the vector of closest proximity. The vector of smallest magnitude that joins the target origin with the missile trajectory is sought. The magnitude of the joining vector is ##EQU18## a quartic function of time through matrix M = M(t). Distance D(t) is minimum at a time t = t min  satisfying ##EQU19## a cubic equation. Determing t min  by means of Newton&#39;s method, the minimum miss distance vector is given by 
     
         MMV = M(t.sub.min) x. 
    
     Of the three trajectory estimation techniques discussed, the Steepest Descent algorithm is the slowest, the One Step Conjugate Gradients Algorithum is between 5 and 10 times faster and the N-step conjugate Gradients Algorithm utilizing 100 steps is approximately 10 times faster than the One Step Conjugate Gradients method. While it is clear, then, that a 100 Step Conjugate Gradients Algorithm is the most efficient, any one of the three systems may be used to solve the problem in the system of the invention. 
     Related conjugate directions algorithms such as the &#34;Davidon-Fletcher-Powell&#34; may be used with equal success. (See Ortega and Rheinboldt, supra, at p. 248.) In general, any descent or relaxation method may be used. 
     Data process 36 may be any one of commercially available computers such as, for example, Xerox Data System Model Sigma 5, properly programmed. The FORTRAN IV program which follows has been used with a Sigma 5 computer in a simulation of data processor 36 and has been found effective. 
     The FORTRAN IV program contains not only the estimation scheme essential to the proper operation of the system, but provisions, as well, for evaluating its performance by means of computational error analysis and generating appropriate statistics. As given, the program consists of several distinct parts: namely, 
     Main : the controlling program that calls primary subroutines RTDATA, SD, MISVEC, and PLOTi. 
     Rtdata : the subroutine that reads in necessary program control parameters and system specifications as well as detected data. In addition, this subroutine perturbs the given data in accordance to an error process for the purpose of evaluating the performance of the estimation procedure with data corrupted by noise. The latter feature is, of course, not essential to the operation of the system. 
     Sd : with detection information available through RTDATA, this subroutine exercises the conjugate gradients estimation process. The end result is a trajectory description in terms of vector acceleration, a, vector velocity, v, and vector displacement, s. 
     Misvec : with a trajectory specified according to SD, this subroutine computes the minimum miss vector, the vector that connects the target origin to the projectile at the time of closest proximity, given in target coordinates. In addition, this subroutine generates appropriate statistics useful in evaluating the performance of the estimation process using noisy detection data. 
     Ploti : This subroutine generates a histogram of vector magnitude errors. It is not essential to the operation of the system. 
     Root : this secondary subroutine serves SD as well as MISVEC in computing roots of cubic equations involved in each. 
     Included also is a listing of 156 data cards necessary for the present program to work as a simulation program, given that appropriate control cards are used. The detected data given has been generated computationally. In actual operation such data will be provided by the system of the invention. A program illustrative of the present invention is set out herein below. ##SPC1##10/924 ##SPC2## 
     While the foregoing description of the preferred embodiment of the invention discloses a scenario in which the miss distance of a missile with respect to a moving airborne target is measured, it will be apparent to one skilled in the art that there may be other applications for the invention. Since the system, as described, computes a trajectory with respect to the &#34;target&#34;, it is of no consequence to the invention if the &#34;target&#34; is not moving. The components of the system of the invention, as herein described as mounted on a &#34;target&#34;, could as well be mounted on a ground or water based mobile vehicle or on such a vehicle in a fixed location or at a fixed (nonmobile) ground based station. The system may be used to accurately record the trajectory of any moving vehicle as well as the missile heretofore described. 
     Various other modifications and changes may be made to the present invention from the principles of the invention described above without departing from the spirit and scope thereof, as encompassed in the accompanying claims.