## 1.

## Introduction

Detection and tracking of fast moving incoming targets are important for a variety of military applications such as active vehicle protection (AVP) systems and close-in weapon systems mounted on modern warships. AVP systems are designed to protect armored military vehicles from incoming threats with reactive weapons for counteractions, so they generally consist of a radar system to track and activate weapons control against incoming high velocity threats.^{1} These threats range from shoulder-launched rocket-propelled grenades to more lethal kinetic energy penetrators. A radar is an active sensor, i.e., it transmits a signal that the target reflects and the reflected signal is collected by the radar’s aperture. This active nature of the radar, while making it an effective sensor, also leaves the vehicle vulnerable to detection and in the worst case, the radar system is jammed by enemy electronic countermeasure systems.^{2} An AVP system with a passive sensor, which measures natural emissions of targets, can be designed to avoid detection and jamming. Considerable research and development have been done to develop algorithms for tracking and position estimation based on radar with laser illumination or infrared (IR) sensors.^{3}^{,}^{4} However, standalone use of multiple IR sensors not only proves to be useful in target tracking and estimation^{5}6.7.^{–}^{8} but also has an added benefit of stealth due to the passive nature.

The bearing only estimation and tracking-based approach involving two or more passive sensors employs nonlinear filter-based approaches like the extended Kalman filter (EKF)^{6}^{,}^{9}^{,}^{10} or the unscented Kalman filter.^{11} These methods use angle-only information measured by multiple systems at the same time. These algorithms are generally used for medium-range target tracking where an accuracy of a few meters is considered sufficient.^{6}^{,}^{9}10.^{–}^{11} An approach for a multiple IR sensor-based system is proposed in Ref. 12, where the estimates of the EKF are evaluated by a fuzzy logic algorithm that decides the authority of 3-D angle intersection results. The scenario described in Ref. 12 is for long range target tracking and the root mean square error (RMSE) results are too large for AVP system application. Likewise, Ref. 7 makes use of the ratio of irradiance measured by two sensors to solve the problem that arises when the flightpath of a target is nearly collinear to the baseline. Although the aforementioned technique improves the overall estimation performance and provides a basic structure for position estimation using measured irradiance, it does not fulfill the requirements of an active protection system for close-range threats. This is because with a small baseline separation between sensors, the measured irradiances of the two sensors will be nearly the same, and consequently the ratio between these irradiances will not provide sufficient information to improve the estimation process. While effective with large baseline and sensor-target distances, the technique presented in Ref. 7 is not suitable for use in AVP systems due to short-range high velocity incoming threats and limitations on the short baseline distance between sensors.

Time-to-go (TTG) is defined as the time required for an object to reach a certain location, e.g., the time a rocket takes from any point of interest in space after being launched to impact at the aimed location. The TTG measurement or estimate can be augmented to angle measurements to improve the overall target state estimation performance. There are two major ways to estimate the TTG from image-based IR sensors, i.e., image processing algorithms and the use of measured target intensity. Estimation from imaging sensors is generally based on optical flow^{13}^{,}^{14} or algorithms for physical feature extraction^{15}^{,}^{16} from time-varying imagery including IR imagery.^{17} The second approach uses the rate of change of intensity to estimate the TTG for an incoming target,^{2} also known as amplitude rate tracking. This technique requires basic image processing on an intensity map of the surveillance region and conversion of each probable object to a point source by obtaining the center of mass of the object as well as the mean intensity. Two different approaches to estimate TTG from amplitude rate tracking are presented in Ref. 2, the first assumes prior knowledge of accurate initial position and parameters related to the irradiance of target types, whereas the second approach considers either already known variables or optimized values of these variables obtained via a trial and error approach. These TTG estimation techniques have certain drawbacks; computational complexity makes the approach involving image processing unsuitable for real-time applications, whereas lack of generality is the main problem in Ref. 2.

The first contribution is the calculation of the Cramér–Rao lower bound (CRLB)^{18} on TTG estimate, which then is used as a covariance for the TTG estimate when used as an additional measurement to the EKF for target localization. The second contribution includes the development of a novel technique to extract the TTG estimate from the target irradiance by using a fixed interval iteration scheme (Newton–Raphson’s method). It is also shown that the resulting TTG estimate can be augmented to angle measurements of the two sensors to improve the position and velocity estimates of the incoming projectile using an EKF in real time. In addition, performance comparison results of the proposed technique with an AVP system based on radar only and another one based on data fusion of radar and an IR sensor are provided.

In the next section, position estimation using multiple passive sensors is introduced along with target measurement models, ranging systems, estimation of the TTG from successive irradiance measurements, and calculation of the CRLB for TTG; then automatic triangulation is explained in Sec. 3. Section 4 describes the TTG estimation and the augmented measurement scenario. The effectiveness of the proposed algorithm for both TTG estimation and reduced RMSE in position is demonstrated along with the performance comparison with a radar-based AVP system via a simulation study in Sec. 5, which is followed by concluding remarks in Sec. 6.

## 2.

## Three-Dimensional Position Estimation with Multiple Passive Sensors

This paper considers multiple IR sensors installed on the same platform arranged in pairs such that each pair of sensors has the same spectral characteristics, shares the same field of view and has perfect alignment. Each pair of sensors will register the same target in its field of view either with the platform moving or stationary. Another classical assumption that we make is the synchronization of each sensor pair. The targets for an AVP system range from subsonic small caliber (105 mm) rocket-propelled grenades with high-explosive anti-tank warheads to supersonic kinetic energy projectiles with larger calibers. The proposed IR sensor resolution, discussed in a later section, aims at the detection of a 100 mm caliber target at a range of 1000 m. The state of the platform is assumed to be known with a certain accuracy and the baseline distance between each sensor pair is fixed. The measurements obtained from an IR sensor are azimuth and elevation angles of a target along with the signal irradiance at a known discrete time indexed by $k$. The radar measurement on the other hand includes the target range information along with azimuth and elevation angles.

A nearly constant velocity model is considered for target dynamics and a nonlinear measurement model is employed for target angles for triangulation. The origin of the tracker coordinate system is set as the position of sensor ${S}_{1}$, as shown in Fig. 1. The distance between the two sensors is the baseline ($b$), and the angles (${\u03f5}_{1}$, ${\u03f5}_{2}$) and (${\eta}_{1}$, ${\eta}_{2}$) are the elevation and azimuth angles as observed by ${S}_{1}$ and ${S}_{2}$, respectively.

## 2.1.

### Target Measurement Models

Measurements of the target are obtained in sequence at each sensor. The measurements consist of angle information, ${\mathbf{z}}_{{s}_{k}}^{\theta}$ (azimuth and elevation), along with irradiance, ${\mathit{z}}_{{s}_{k}}^{w}$, from both sensors. Irradiance is the ratio of the intensity of the radiation emitted by the target to the square of the distance separating the target from the sensor system with losses due to atmospheric absorption and scattering accounted for.^{2} The measurement model for azimuth and elevation angles with the relative Cartesian state vector ${\mathbf{X}}_{k}$ can be expressed as

## (1)

$${\mathbf{z}}_{{s}_{k}}^{\theta}={\mathbf{h}}_{k}^{\theta}({\mathbf{X}}_{k})+{\mathbf{w}}_{k}^{\theta}.$$The subscript $s$ denotes the sensor number generating measurement, where ${\mathbf{w}}_{k}^{\theta}$ is the zero-mean, white Gaussian measurement noise with covariance ${\mathbf{R}}^{\theta}$, which is expressed by

## (2)

$${\mathbf{R}}^{\theta}\u2254\mathrm{diag}({\sigma}_{{\eta}_{1}}^{2},{\sigma}_{{\epsilon}_{1}}^{2},{\sigma}_{{\eta}_{2}}^{2},{\sigma}_{{\epsilon}_{2}}^{2})$$## (3)

$${\mathbf{h}}_{k}^{\theta}({X}_{k})={[\begin{array}{cccc}{\eta}_{{1}_{k}}& {\epsilon}_{{1}_{k}}& {\eta}_{{2}_{k}}& {\epsilon}_{{2}_{k}}\end{array}]}^{\prime},$$## (4)

$${\eta}_{{s}_{k}}={\mathrm{tan}}^{-1}\frac{{y}_{k}^{t}-{y}_{k}^{s}}{{x}_{k}^{t}-{x}_{k}^{s}},\phantom{\rule{0ex}{0ex}}{\epsilon}_{{s}_{k}}={\mathrm{tan}}^{-1}\frac{{z}_{k}^{t}-{z}_{k}^{s}}{{({({x}_{k}^{t}-{x}_{k}^{s})}^{2}+{({y}_{k}^{t}-{y}_{k}^{s})}^{2})}^{1/2}}.$$The following relationship was used to relate irradiance, intensity, and the distance of sensor ${S}_{1}$ from the target^{2}

^{19}where $\alpha \u200a({m}^{-1})$ is the extinction coefficient and can be calculated as and the maximum value of $\mathrm{\Delta}b$ does not exceed 0.3 $(\mathrm{\Delta}b\le 0.3)$. The measurement thus obtained can be written as

## (7)

$${\mathit{z}}_{{s}_{k}}^{w}=\frac{{I}_{o}}{{r}_{{s}_{k}}^{2}}{\mathrm{e}}^{-\alpha {r}_{k}}+{\mathit{n}}_{k},$$## (8)

$${\mathit{z}}_{{s}_{k}}^{w}=\frac{{I}_{o}}{{v}_{k}^{2}{({t}_{f}-kT)}^{2}}{\mathrm{e}}^{-\alpha {v}_{k}({t}_{f}-kT)}+{\mathit{n}}_{k},$$Instead of using the measured irradiance ${\mathit{z}}_{{s}_{k}}^{w}$ directly, the measured irradiance is used to calculate the TTG, a procedure that forms the basis of this research, and will be explained in Sec. 2.3.

## 2.2.

### Ranging Systems

A key function of combining angle-only information from two passive sensors is to be able to observe the range of the target. The location of a point is determined by measuring the angles to it from two points separated by a baseline, as shown in Fig. 1. Since this paper deals with the ranging issue, the problem is slightly modified to emphasize the target’s position or the distance from the target to the expected point of impact. The origin of the tracking coordinate system in our case is located at the position of sensor 1, as shown in Fig. 1. The ground range $a$ can be expressed as

## (9)

$$a=\frac{b\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}{\eta}_{2}}{\mathrm{sin}({\eta}_{2}-{\eta}_{1})},$$There are different ways to address this geometric problem; using Eqs. (9) and (10) is one of them.^{7} The position of the target calculated from the baseline length and angle information is used in this paper to initialize the triangulation process. In order to estimate the target position, angles from both sensors separated by a fixed baseline are augmented together and the estimation process is accomplished using the EKF.^{7}^{,}^{20}^{,}^{21}

Although the process of initialization using Eqs. (9) and (10) seems straightforward and very accurate due to noisy angle measurements from the two passive sensors in Fig. 1 and prediction errors, the sensor-to-target vectors may not intersect at the target position.^{2} Using the known noise variance in angle measurement, a region of interest is defined (indicated by the shaded region in Fig. 2). Based on this region, the statistical parameters of Gaussian noise affecting the true target range are defined.

The parameters ${\eta}^{+}$, ${\eta}^{-}$, ${\epsilon}^{+}$, and ${\epsilon}^{-}$ are the three sigma bounds on the measurement noise for the azimuth and elevation angles, as in Eq. (2). The parameter $\delta $ varies in each scan, which implies the variance of the noise being added in the target range will change at every scan. The variance of this statistical parameter for range is defined as

Bearing in mind that the cross-range noise variance is considerably smaller compared to the down-range noise variance, the latter will be the dominating source of RMSE in target range estimation.

## 2.3.

### Time-To-Go Estimation and Cramér–Rao Lower Bound

To exploit the full potential of the passive IR sensors, the measured irradiance from these sensors is used to estimate the TTG. An initial guessed value of the parameter is needed for initializing the recursive estimation algorithm based on the Newton–Raphson method for parameter estimation. The irradiance measurement in Eq. (8) is obtained from two consecutive scans and then the ratio between these measurements is used to get a sufficiently close value for initialization of the estimator for the TTG. The measurement Eq. (8) for $k=\mathrm{0,1}$ can be rewritten as

## (12)

$${\mathit{z}}_{{s}_{0}}^{w}=\frac{{A}_{v}{\mathrm{e}}^{-\alpha {v}_{k}{t}_{f}}}{{t}_{f}^{2}}+{n}_{0},$$## (13)

$${\mathit{z}}_{{s}_{1}}^{w}=\frac{{A}_{v}{\mathrm{e}}^{-\alpha {v}_{k}({t}_{f}-T)}}{{({t}_{f}-T)}^{2}}+{n}_{1},$$^{22}After analysis and carefully carried out simulations for convergence, it was concluded that in the calculation of the inital value for the TTG, the irradiance measurement noise can be ignored for fast convergence by the Newton–Raphson algorithm. Thus, the ratio of the irradiances obtained from the same sensor at different scans, after ignoring the measurement noise irradiance can be expressed as

## (14)

$$\gamma =\frac{{\mathit{z}}_{{s}_{1}}^{w}}{{\mathit{z}}_{{s}_{0}}^{w}}=\frac{{t}_{f}^{2}{\mathrm{e}}^{\alpha {v}_{k}T}}{{({t}_{f}-T)}^{2}}.$$It was also concluded that for the system range variation in the engagement environment under consideration, i.e., less than 1000 m, the effects of signal (irradiance) attenuation due to atmosphere are negligible and can be ignored. This assumption is made for rest of the paper, i.e., the term, ${\mathrm{e}}^{\alpha {v}_{k}T}$ will be approximated by one.^{8} So with this assumption, Eq. (14) becomes

After calculating the value of ${t}_{f}$, Eq. (15) is used to obtain the initial value for the other unknown parameter, ${A}_{v}$, such that

## (16)

$${A}_{v}={\mathit{z}}_{{s}_{0}}^{w}{\left(\frac{\sqrt{\gamma}T}{\sqrt{\gamma}-1}\right)}^{2}.$$Suppose an unbiased estimate of the parameter set $\theta $ composed of ${t}_{f}$ and the parameter ${A}_{v}$ can be generated from the simplified version of Eq. (8). Then the estimation error variance will be bounded by the inverse of the Fisher information matrix, i.e., $E[{(\widehat{\theta}-\theta )}^{2}]\ge {I}^{-1}(\theta )$. The information matrix can be computed by taking the expected value of the Hessian of the log-likelihood function as follows:

## (17)

$$\mathbf{I}(\theta )=E\{{[{\nabla}_{\theta}\text{\hspace{0.17em}}\mathrm{log}\text{\hspace{0.17em}}p({z}_{{s}_{k}}^{w}|\theta )]}^{T}[{\nabla}_{\theta}\text{\hspace{0.17em}}\mathrm{log}\text{\hspace{0.17em}}p({z}_{{s}_{k}}^{w}|\theta )]\}.$$The likelihood function that will be used to calculate the information matrix can be written as

## (18)

$$p({z}_{{s}_{k}}^{w}|\theta )=\mathcal{N}({z}_{{s}_{k}}^{w};{\mathit{h}}_{k}^{w},{\sigma}_{w}^{2}),$$## (19)

$$p({\mathbf{Z}}_{s}^{w}|\theta )=\prod _{k=0}^{n-1}\mathcal{N}({z}_{{s}_{k}}^{w};{\mathit{h}}_{k}^{w},{\sigma}_{w}^{2}),\phantom{\rule{0ex}{0ex}}=\frac{1}{{(2\pi )}^{n/2}{\sigma}_{w}^{n}}{\mathrm{e}}^{-\frac{1}{2{\sigma}_{w}^{2}}\sum _{k=0}^{n-1}{({z}_{{s}_{k}}^{w}-{h}_{k}^{w})}^{2}}.$$By the assumption of negligible atmospheric attenuation, the measurement function becomes

## (20)

$${h}_{k}^{w}=\frac{{I}_{o}}{{v}^{2}}\frac{1}{{({t}_{f}-kT)}^{2}}=\frac{{A}_{v}}{{({t}_{f}-kT)}^{2}},$$## (21)

$$P=\mathrm{log}\text{\hspace{0.17em}}p({\mathbf{Z}}_{s}^{w}|\theta ),\phantom{\rule{0ex}{0ex}}=-\mathrm{log}({(2\pi )}^{n/2}{\sigma}_{w}^{n})-\frac{1}{2{\sigma}_{w}^{2}}\sum _{k=0}^{n-1}{({z}_{{s}_{k}}^{w}-{h}_{k}^{w})}^{2},$$## (22)

$$\frac{\partial P}{\partial {t}_{f}}=\frac{1}{2{\sigma}_{w}^{2}}\sum _{k=0}^{n-1}2({z}_{{s}_{k}}^{w}-{h}_{k}^{w})\left(\frac{\partial {h}_{k}^{w}}{\partial {t}_{f}}\right),$$## (23)

$$\frac{\partial P}{\partial {t}_{f}}=\frac{-2{A}_{v}}{{\sigma}_{w}^{2}}\sum _{k=0}^{n-1}\frac{({z}_{{s}_{k}}^{w}-{h}_{k}^{w})}{{({t}_{f}-kT)}^{3}},$$## (24)

$$\frac{\partial P}{\partial {A}_{v}}=\frac{1}{{\sigma}_{w}^{2}}\sum _{k=0}^{n-1}\frac{({z}_{{s}_{k}}^{w}-{h}_{k}^{w})}{{({t}_{f}-kT)}^{2}}.$$By combining Eqs. (23) and (24), the gradient of the log-likelihood function with respect to the vector $\theta $ can now be written as

## (25)

$${\nabla}_{\theta}\text{\hspace{0.17em}}\mathrm{log}\text{\hspace{0.17em}}p({\mathbf{Z}}_{s}^{w}|\theta )={\left[\begin{array}{c}\frac{-2{A}_{v}}{{\sigma}_{w}^{2}}\sum _{k=0}^{n-1}\frac{({z}_{{s}_{k}}^{w}-{h}_{k}^{w})}{{({t}_{f}-kT)}^{3}}\\ \u200a\u200a\frac{1}{{\sigma}_{w}^{2}}\sum _{k=0}^{n-1}\frac{({z}_{{s}_{k}}^{w}-{h}_{k}^{w})}{{({t}_{f}-kT)}^{2}}\end{array}\right]}^{\prime}.$$From Eqs. (17) and (25), a $2\times 2$ matrix is obtained such that the first and second diagonal entries of the inverse matrix are the CRLBs for TTG and ${A}_{v}$. Then the information matrix can be written as

## (26)

$$\mathbf{I}(\theta )=\left(\begin{array}{cc}E\{{D}_{11}\}& E\{{D}_{12}\}\\ E\{{D}_{21}\}& E\{{D}_{22}\}\end{array}\right),$$## (27)

$$E\{{D}_{11}\}=\frac{4{A}_{v}^{2}}{{\sigma}_{w}^{2}}\sum _{k=0}^{n-1}\frac{1}{{({t}_{f}-kT)}^{6}},$$## (28)

$$E\{{D}_{12}\}=E\{{D}_{21}\}=\frac{-2{A}_{v}}{{\sigma}_{w}^{2}}\sum _{k=0}^{n-1}\frac{1}{{({t}_{f}-kT)}^{5}},$$Finally, the CRLB is calculated using Eqs. (27)–(29) in Eq. (26) and then inverting the Fisher information matrix. Then

where the first and second diagonal entries of Eq. (30) represent the CRLB of the desired parameters and ${\mathbf{P}}_{0|0}$ is the initial covariance matrix.The CRLB computed here is used to determine whether the TTG estimate is useful in reducing the position RMSE at a particular signal-to-noise ratio (SNR) of the target irradiance and background noise; the second purpose is to use the CRLB as the measurement noise of the TTG estimate when used as an augmented measurement to the angle measurements of the EKF for target localization with the IR sensors.

## 3.

## Ranging Based on Automatic Triangulation

Due to the stochastic nature of the measurements from the two sensors, the EKF is used for target state estimation.^{18} The target’s trajectory state propagates in the discrete domain by

The state transition matrix $\mathbf{F}$ in Eq. (31) is defined as

## (33)

$$F=\left(\begin{array}{cc}{\mathbf{I}}_{3}& T{\mathbf{I}}_{3}\\ {\mathbf{O}}_{3}& {\mathbf{I}}_{3}\end{array}\right),{\mathbf{I}}_{3}=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right),$$## (34)

$$\mathbf{Q}=q.\left(\begin{array}{cc}\frac{{T}^{4}}{4}{I}_{3}& \frac{{T}^{3}}{2}{I}_{3}\\ \frac{{T}^{3}}{2}{I}_{3}& {T}^{2}{I}_{3}\end{array}\right),$$^{2}

^{,}

^{18}can be expressed as

## (36)

$${\widehat{\mathbf{P}}}_{k/k-1}={\mathbf{F}\widehat{\mathbf{P}}}_{k-1/k-1}{\mathbf{F}}^{T}+\mathbf{Q},$$## (37)

$${\mathbf{S}}_{k}={\mathbf{H}}_{k}^{\theta}{\widehat{\mathbf{P}}}_{k/k-1}{({\mathbf{H}}_{k}^{\theta})}^{T}+{\mathbf{R}}^{\theta},$$## (38)

$${\mathbf{K}}_{k}={\widehat{\mathbf{P}}}_{k/k-1}{({\mathbf{H}}_{k}^{\theta})}^{T}{\mathbf{S}}_{k}^{-1},$$## (39)

$${\widehat{\mathbf{X}}}_{k/k}={\widehat{\mathbf{X}}}_{k/k-1}+{\mathbf{K}}_{k}[{\mathbf{z}}_{{(\mathrm{1,2})}_{k}}^{\theta}-{\mathbf{h}}_{k}^{\theta}({\widehat{\mathbf{X}}}_{k/k-1})],$$## (40)

$${\widehat{\mathbf{P}}}_{k/k}=(\mathbf{I}-{\mathbf{K}}_{k}{\mathbf{H}}_{k}^{\theta}){\widehat{\mathbf{P}}}_{k/k-1},$$## (41)

$${\mathbf{H}}_{k}^{\theta}=\left(\begin{array}{cccccc}\frac{\partial {\eta}_{{1}_{k}}}{\partial x}& \frac{\partial {\eta}_{{1}_{k}}}{\partial y}& 0& 0& 0& 0\\ \frac{\partial {\epsilon}_{{1}_{k}}}{\partial x}& \frac{\partial {\epsilon}_{{1}_{k}}}{\partial y}& \frac{\partial {\epsilon}_{{1}_{k}}}{\partial z}& 0& 0& 0\\ \frac{\partial {\eta}_{{2}_{k}}}{\partial x}& \frac{\partial {\eta}_{{2}_{k}}}{\partial y}& 0& 0& 0& 0\\ \frac{\partial {\epsilon}_{{2}_{k}}}{\partial x}& \frac{\partial {\epsilon}_{{2}_{k}}}{\partial y}& \frac{\partial {\epsilon}_{{2}_{k}}}{\partial z}& 0& 0& 0\end{array}\right).$$Each element in Eq. (41) can be expressed using partial derivatives of Eq. (4) in the following generalized form (for both sensors), where the time index has been removed for simplicity and the superscript denotes the sensor number generating the measurement

The Cartesian coordinates are not known and are written in a simplified manner here; the actual elements of the Jacobian matrix are calculated using the state vector obtained after prediction. The filter covariance matrix is initialized based on Eq. (11). The quality of the estimate depends on the sampling time $T$ along with the assumed measurement process noise covariance matrix ${\mathbf{R}}^{\theta}$ and process noise variance $q$.

## 4.

## Ranging Based on Triangulation and Time-To-Go Estimation

The main emphasis of this work is to show how to estimate TTG in an iterative manner based on intensity information from a passive IR sensor at short ranges with a high degree of accuracy. The case where estimation is performed using triangulation with an EKF results in a performance that needs improvement, both with respect to response time and the RMSE for AVP. With the high speed projectiles and shorter ranges under consideration, time is of the essence. Not only does the RMSE in position need to be decreased, but it has to be done quickly.

A well-known fixed interval iteration algorithm, Newton–Raphson, is used for estimating the TTG. The initialization process for the case of triangulation at greater ranges causes the EKF to diverge, which is because the shaded area representing possible target location as indicated in Fig. 2 is very large, thus it produces unacceptable values for initialization of the EKF. With the noise statistics for angle measurements being considered in this research, filter divergence occurs at target ranges greater than 600 m. This discussion is based on a baseline separation of 3 m between the two sensors, the resolution of each sensor is $1280\times 1024$ pixels and the target velocity is $1000\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}/\mathrm{s}$. The irradiance measurements on the other hand can still be used to estimate the TTG beyond this range. The triangulation process is not started until the range of the projectile is less than 600 m with the uncertainty in range accounted for. If the projectile is fired from a range greater than this bound then all the measurements from the detection up to this point in time are used to estimate the TTG. During this period, the sampling time for the irradiance measurement is kept intentionally large. This strategy is adopted so that the initialization part, i.e., Eq. (15), gives a close enough initial value of ${t}_{f}$ because of the higher levels of noise at these ranges. Violating this strategy results in an increased number of the Newton–Raphson iterations, thus exceeding the predefined limit. Otherwise the estimator fails to converge. As the target moves closer, more samples are drawn out and used to estimate the TTG because of the improved target irradiance-to-background noise ratio. Specifically, when the range is less than 500 m, the TTG is estimated based on target irradiance measured over every 100 m. The EKF with the augmented measurement vector starts working as soon as the first TTG estimate is available.

The measurement augmentation takes place only after a TTG estimate is available. Since target tracking algorithms also need some time to confirm a track, the fast sampling rate in this case provides sufficient measurements, making the track confirmation process quicker. The algorithm here switches between two filters, i.e., an EKF for position estimation using Eqs. (9) and (10) for initialization and an EKF with an extended measurement vector, the additional measurement being the estimated TTG. The complete flow of this algorithm is given in Fig. 3.

## 4.1.

### Time-To-Go Estimation

Using Eq. (20) as the measurement generating function and differentiating with respect to $\theta $, which is a column vector containing the parameters to be estimated as defined earlier, one has

## (42)

$$\frac{\partial {h}_{k}^{w}}{\partial {\theta}_{N}}=(\frac{-2{A}_{v}}{{({t}_{f}-kT)}^{3}},\frac{1}{{({t}_{f}-kT)}^{2}}).$$Table 1 presents a pseudocode of the Newton–Raphson method for TTG estimation. The algorithm is initialized using Eqs. (15) and (16). After careful analysis and simulation, the number of iterations was limited to 10, since the algorithm converges in less than 5 iterations.

## Table 1

Newton-Raphson for time-to-go estimation.

Input: Irradiance measurement ${\mathit{z}}_{{\mathit{s}}_{k}}^{\mathit{w}}$ |

if range less than bound then |

Gather minimum required input samples. |

if required samples criteria are met |

calculate $\theta $ for start point using (15) and (16) |

for number of iterations do |

Calculate: |

$M=\sum _{k=0}^{n-1}{(\partial {h}_{\mathrm{k}}^{\mathrm{w}}/\partial {\theta}_{N})}^{T}{\sigma}_{w}^{-2}(\partial {h}_{\mathrm{k}}^{\mathrm{w}}/\partial {\theta}_{N})$ |

$N=\sum _{k=0}^{n-1}{(\frac{\partial {h}_{\mathrm{k}}^{\mathrm{w}}}{\partial {\theta}_{N}})}^{T}{\sigma}_{w}^{-2}({z}_{{\mathrm{s}}_{\mathrm{k}}}^{\mathrm{w}}-{h}_{\mathrm{k}}^{w})$ |

Update: |

${\widehat{\theta}}_{n}^{l+1}={\widehat{\theta}}_{n}^{l}+{M}^{-1}N$ |

end for |

Output:Estimated parameters, ${t}_{f}$ and $A$, |

end if |

else |

Gather minimum required input samples at low sampling rate |

end if |

## 4.2.

### Filter Structure with Extended Measurement Vector

The EKF is employed again, with the state vector consisting of position and velocity as in the previous case but with an extended measurement vector ${\mathbf{Z}}_{{s}_{k}}$ and hence an extended measurement matrix ${\mathbf{H}}_{k}$. The new measurement vector, which consists of five elements, i.e., target azimuth and elevation from two sensors along with TTG $({\tau}_{k})$, can be expressed as

## (43)

$${\mathbf{Z}}_{{s}_{k}}={[\begin{array}{cccc}{\eta}_{{1}_{k}}& {\epsilon}_{{1}_{k}}& {\eta}_{{2}_{k}}& {\epsilon}_{{2}_{k}}\end{array}]}^{T}.$$The measurement matrix will have an additional row of partial derivatives of $({\tau}_{k})$ with respect to the state vector, which can be written as

## (44)

$${\mathbf{H}}_{k}^{\tau}=\left[\begin{array}{cccccc}\frac{\partial {\tau}_{k}}{\partial x}& \frac{\partial {\tau}_{k}}{\partial y}& \frac{\partial {\tau}_{k}}{\partial z}& \frac{\partial {\tau}_{k}}{\partial \dot{x}}& \frac{\partial {\tau}_{k}}{\partial \dot{y}}& \frac{\partial {\tau}_{k}}{\partial \dot{z}}\end{array}\right],$$## (46)

$${\mathbf{H}}_{k}=\left[\begin{array}{c}{\mathbf{H}}_{k}^{\theta}\\ {\mathbf{H}}_{k}^{\tau}\end{array}\right].$$Similarly, the corresponding covariance matrix needs to be extended with the last appended diagonal entry being the variance for TTG. The TTG estimate obtained from Table 1 is used as an additional input to EKF and is considered to be noisy, where the additive noise follows a normal distribution with zero mean and its variance equals the value obtained from Eq. (27). The variance changes for TTG at every point in time when the extended measurement vector is used, because, with the target getting closer and closer to the sensors, the irradiance measurement becomes less noisy.

Except for Eq. (46) and the extended measurement noise matrix, the rest of the structure is the same as the classical EKF defined by Eqs. (35)–(40). The estimation process continues until it reaches a minimum threshold distance for the target to be engaged.

## 5.

## Simulation Study

The estimation performance for both TTG and target localization is provided in terms of RMSE. The simulation scenario is depicted in Fig. 1. Two passive IR sensors are placed along the positive $x$-axis with the first sensor positioned at [0, 0, 2 m] and the second one at [$b$, 0, 2 m]. The baseline distance $b$ is expressed in meters. The target is moving almost parallel to the $y$-axis following a uniform motion model with velocity $v=1000\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}/\mathrm{s}$. The initial position of the target is determined by the vector [1.5 m, 1000 m,2 m], with the impact point centered between the two sensors.

There are two criteria for optimal baseline length selection. The first one is the variance of target range estimate and the second one is the availability of space on the vehicle for installation of IR sensors. The first criterion affects the overall estimation performance, as a larger variance means a larger RMSE in position, which cannot be used for an AVP system. The second criterion, i.e., availability of space on the vehicle limits the baseline length physically to a maximum of $b=3\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}$. As illustrated in Fig. 4, the variance of the range with a baseline length of $b=1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}$ calculated based on Eq. (11) is substantially large and the filter fails to converge when estimating the target position within the limited time available. With these results, we chose the maximum allowable baseline distance $b=3\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}$.

The sampling time for the triangulation scheme was set to 0.01 s, whereas in the case of TTG estimation it was varied from 0.1 s before the first target state estimation to 0.01 s thereafter. The sensors considered for this scenario have a resolution of $1280\times 1024$ pixels, with an aspect ratio of 5:4 and a field of view (FOV) of 7.3 deg along azimuth, and 5.9 deg elevation. However, this requirement can be relaxed based on the target detection range or by narrowing the FOV. With the current high-resolution sensors,^{23}^{,}^{24} it is reasonable to use sensors with the resolution and sampling time considered for this simulation scenario. The maximum angle deviation for each sensor is kept at one pixel, which corresponds to ${\sigma}_{\eta}=5.7\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mrad}$ and ${\sigma}_{\epsilon}=5.7\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mrad}$.

The simulation in case of the AVP system based on IR sensors only was conducted for different values of SNR between the target irradiance given by Eq. (5) and background noise, ${n}_{k}$. While the CRLB suggested an SNR of up to 15 dB will result in improved position estimation after measurement augmentation, filter divergence limits the SNR to 25 dB. This was due to the unsuitable initial points obtained from Eqs. (15) to (16). The CRLBs for different values of SNR at TTG from $1-0.1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{s}$ are given by Figs. 5 and 6.

We also compare the proposed approach for the AVP system with already existing approaches based on radar and radar-IR fusion. The sensor location, in case of the AVP system based on radar, is set to be at [0, 0, 2 m]. The uncertainty levels in the target measurements of the radar are defined by the diagonal elements of the measurement covariance matrix, ${\mathbf{R}}^{R}$, and are quantified by standard deviations, ${\sigma}_{r}=1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}$, ${\sigma}_{\epsilon}=26.4\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mrad}$, and ${\sigma}_{\eta}=26.4\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mrad}$. The target has the same velocity while the initial position vector is determined by [0, 1000 m,2 m]. The sampling time is equal to that of the IR sensors, i.e, $T=0.01\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{s}$. A Kalman filter is used to estimate the position and velocity of the target by utilizing the unbiased converted measurements in Cartesian coordinates as described in Refs. 2 and 25. In the next case where the radar measurements and the measurements of an IR sensor are fused, the radar location is not changed, but the IR sensor location is given by the vector [0.5 m, 0, 2 m]. The target position is also kept the same as in the case of the radar based system. The IR sensor used here shares the same specifications and statistical parameters as in the case of an IR only AVP system. The predicted target impact point in this case is the radar position. It is also assumed that both the sensor systems are perfectly synchronized. The Kalman filter is used to update the state estimate with the converted radar measurements whereas an EKF is used to update the IR information by utilizing the updated Kalman filter state estimates as the predicted state estimates.

In order to demonstrate performance improvement in the ranging system and the efficiency of TTG estimation, Monte Carlo simulation results are presented. Each simulation experiment of triangulation with TTG, radar only and radar-IR fusion consists of 1000 runs.

The following figures are presented to illustrate and compare the performance of the algorithms. Figure 7(a) shows the estimated TTG based on the Newton–Raphson algorithm where each estimate is obtained as a result of 10 iterations. The TTG was estimated in two parts as illustrated; the first estimate was made using measurements from 1000 to 500 m with a sampling rate of 0.1 s, the following batch was estimated after every 10 measurements with a sampling rate of 0.01 s. Figure 7(b) illustrates the quick convergence of the Newton–Raphson algorithm at 100 m for different values of SNR. The estimator here yields near optimal results as compared to the CRLB.

Since the TTG estimate is very accurate, it is bound to improve performance of the EKF. The RMSEs in position and velocity estimates for the systems in comparison are depicted in Figs. 8 and 9, respectively.

The system based on IR sensors only gives a better performance as compared to the system based on the radar and radar-IR fusion approach because of the more accurate angle measurements generated by the IR sensors. Even though the radar-IR fusion approach seems to give the best performance regarding position estimation up to $\mathrm{TTG}=0.4\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{s}$, the passive sensor systems show a better performance afterward. Both systems based on radar only and radar-IR fusion show the position estimates with an RMSE bigger than 1.5 m at the target range of 100 m, which is not accurate enough to engage the target with the counteracting gun system. The bad performance of radar-based systems is due to the dominating factor of the radar range covariance, ${\sigma}_{r}$, and inherent high measurement noise associated with the angle measurements. The radar-IR approach gives an improved performance initially but fails to compete with the passive systems after $\mathrm{TTG}=0.25\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{s}$. A similar trend can also be seen in the RMSE in the velocity estimatation. Therefore, the performance of passive systems with or without TTG augmentation is superior to that of the radar systems. In case of the two IR sensors based approaches, the improvement is clear when the RMSE obtained by the two techniques is compared. The RMSE in the position for triangulation at 100 m is approximately equal to 0.25 m, whereas with the TTG information augmentation, this falls to 0.03 m, which accounts for approximately a 10-fold improvement in results. Reactive weapons to destroy incoming threats are usually deployed within 30 to 50 m from the AVP system, thus a decrease in the down-range RMSE is expected.

To simulate one scan time of 10 ms, the average CPU execution time of the proposed algorithm is 1.7 ms, whereas 1.85 ms is required for the case of radar-IR fusion, which is slightly longer than the proposed algorithm. In the radar only case, the execution time is 0.9 ms. The execution time is calculated for the MATLAB 8.3 platform (3.1 GHz, Intel, Core i5 CPU).

## 6.

## Conclusions

This paper considers the problem of an AVP system based on two passive IR sensors separated by a short baseline distance.

Based on the principle of triangulation using a classical EKF, a new algorithm is devised. The proposed method makes use of a fixed-interval iterative method for estimating TTG and ultimately uses it as an additional TTG measurement for the EKF to improve range estimation. The result is a simple, efficient algorithm with superior estimation performance to the conventional triangulation approach and approaches based on radar or radar-IR fusion. This also results in improved interception at extended range due to a reduction in position RMSE with the capability of the system to prioritize closely fired projectiles with different velocities. Furthermore, the system can serve as a backup in case the systems using radar are jammed or damaged. The simulations show that the accuracy, quick convergence, and simplicity of the proposed algorithm make it suitable for real tactical situations.

The performance of an active protection system based on multiple, closely spaced passive IR sensors can be significantly enhanced with the proposed TTG estimation algorithm against short-range high velocity incoming threats.

## Acknowledgments

This work was supported by the Defense Acquisition Program Administration and the Agency for Defense Development under the contract UD140081CD.

## References

## Biography

**Ihsan Ullah** received his MSc degree in electrical engineering from National University of Science and Technology, Islamabad, Pakistan, in 2010. He was a member of the academic staff at COMSATS Institute of Information Technology, Pakistan, between 2011 and 2012. He is currently working toward his PhD at Hanyang University, Republic of Korea. His research interests include target tracking, estimation, and sensor information fusion.

**Taek Lyul Song** received his BSc degree in nuclear engineering from Seoul National University, Republic of Korea, in 1974 and his MSc and PhD degrees in aerospace engineering from the University of Texas at Austin, USA, in 1981 and 1983, respectively. From 1974 to 1994, he was with the Agency for Defense, Republic of Korea. He has been a professor at the Electronic Systems Engineering Department, Hanyang University, Republic of Korea, since 1995. His research interests include target tracking, guidance, navigation, and control.

**Thia Kirubarajan** received his BA and MA degrees in electrical and information engineering from Cambridge University, England, UK, in 1991 and 1993, and the MS and PhD degrees in electrical engineering from the University of Connecticut, Storrs, USA, in 1995 and 1998. Currently, he is a professor in the Electrical and Computer Engineering Department at McMaster University, Ontario, Canada. His research interests are in estimation, target tracking, multisource information fusion, signal detection, and fault diagnosis.