Active vehicle protection using angle and time-to-go information from high-resolution infrared sensors

Ihsan Ullah; Taek Lyul Song; Thia Kirubarajan

doi:10.1117/1.OE.54.5.053110

13 May 2015 Active vehicle protection using angle and time-to-go information from high-resolution infrared sensors

Ihsan Ullah, Taek Lyul Song, Thia Kirubarajan

Author Affiliations +

Optical Engineering, Vol. 54, Issue 5, 053110 (May 2015). https://doi.org/10.1117/1.OE.54.5.053110

Abstract

An algorithm for an active protection system with two or more passive infrared sensors based on the principle of triangulation and an extended measurement vector consisting of time-to-go (TTG) estimation is proposed to address high-speed short-range threats. The proposed algorithm exploits the full potential of passive sensors enabling the use of target irradiance in a manner to estimate the TTG along with angle information. The approach here uses target irradiance measured at successive scans to initialize a fixed interval iteration scheme for estimating the TTG. The Newton’s method is used for estimating the TTG, the estimate is then added to the measurement vector and position estimation is performed using an extended Kalman filter. Computer simulations demonstrate the resulting improvement in estimation and feasibility for real-time application of the proposed algorithm.

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.¹ 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.² 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.³^,⁴ However, standalone use of multiple IR sensors not only proves to be useful in target tracking and estimation⁵^–⁸ 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)⁶^,⁹^,¹⁰ or the unscented Kalman filter.¹¹ 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.⁶^,⁹^–¹¹ 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¹³^,¹⁴ or algorithms for physical feature extraction¹⁵^,¹⁶ from time-varying imagery including IR imagery.¹⁷ The second approach uses the rate of change of intensity to estimate the TTG for an incoming target,² 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)¹⁸ 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 ( $ϵ_{1}$ , $ϵ_{2}$ ) and ( $η_{1}$ , $η_{2}$ ) are the elevation and azimuth angles as observed by $S_{1}$ and $S_{2}$ , respectively.

Fig. 1

Tracker coordinate system, azimuth angle. $η \in [0,2 π]$ and elevation $ϵ \in [- π / 2, π / 2]$ .

2.1.

Target Measurement Models

Measurements of the target are obtained in sequence at each sensor. The measurements consist of angle information, $z_{s_{k}}^{θ}$ (azimuth and elevation), along with irradiance, $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.² The measurement model for azimuth and elevation angles with the relative Cartesian state vector $X_{k}$ can be expressed as

Eq. (1)

z_{s_{k}}^{θ} = h_{k}^{θ} (X_{k}) + w_{k}^{θ} .

The subscript $s$ denotes the sensor number generating measurement, where $w_{k}^{θ}$ is the zero-mean, white Gaussian measurement noise with covariance $R^{θ}$ , which is expressed by

Eq. (2)

R^{θ} ≔ diag (σ_{η_{1}}^{2}, σ_{ε_{1}}^{2}, σ_{η_{2}}^{2}, σ_{ε_{2}}^{2})

and

Eq. (3)

h_{k}^{θ} (X_{k}) = {[\begin{matrix} η_{1_{k}} & ε_{1_{k}} & η_{2_{k}} & ε_{2_{k}} \end{matrix}]}^{'},

where the bearings and elevation angles are defined as

Eq. (4)

η_{s_{k}} = \tan^{- 1} \frac{y_{k}^{t} - y_{k}^{s}}{x_{k}^{t} - x_{k}^{s}}, ε_{s_{k}} = \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²

Eq. (5)

W = \frac{I_{o}}{r^{2}} e^{- α r},

where

W

, the irradiance, is calculated when the target emits radiation of intensity

I_{o}

at distance

r

. The emitted radiation of intensity,

I_{o}

, mentioned here or in any of the following sections refers to the intensity values associated with the point target on a single frame of the IR image. The exponential term can be considered a sufficient approximation for atmospheric transmittance,¹⁹ where

α (m^{- 1})

is the extinction coefficient and can be calculated as

Eq. (6)

α = \frac{Δ b}{1000} \ln (2),

and the maximum value of

Δ b

does not exceed 0.3

(Δ b \leq 0.3)

. The measurement thus obtained can be written as

Eq. (7)

z_{s_{k}}^{w} = \frac{I_{o}}{r_{s_{k}}^{2}} e^{- α r_{k}} + n_{k},

with

n_{k}

describing the background noise irradiance as in Ref. 7. This noise term appears because the target is viewed against a background which also emits IR energy, e.g., clouds or trees. The denominator can further be expressed such that Eq. (7) becomes

Eq. (8)

z_{s_{k}}^{w} = \frac{I_{o}}{v_{k}^{2} {(t_{f} - k T)}^{2}} e^{- α v_{k} (t_{f} - k T)} + n_{k},

where

v_{k}

represents the target velocity at discrete time

k

,

t_{f}

is the final time, and

T

is the sensor sampling time.

Instead of using the measured irradiance $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

Eq. (9)

a = \frac{b \sin η_{2}}{\sin (η_{2} - η_{1})},

this then can be used to calculate the actual range of the target from the origin as

Eq. (10)

d_{1} = \frac{a}{\cos (ε_{1})} .

There are different ways to address this geometric problem; using Eqs. (9) and (10) is one of them.⁷ 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.⁷^,²⁰^,²¹

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.² 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.

Fig. 2

Geometric relationship between two passive sensors and statistical parameters for Gaussian noise affecting the target range.

The parameters $η^{+}$ , $η^{-}$ , $ε^{+}$ , and $ε^{-}$ are the three sigma bounds on the measurement noise for the azimuth and elevation angles, as in Eq. (2). The parameter $δ$ 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

Eq. (11)

σ_{R_{k}}^{2} = {(\frac{δ_{k}}{3})}^{2} .

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 = 0,1$ can be rewritten as

Eq. (12)

z_{s_{0}}^{w} = \frac{A_{v} e^{- α v_{k} t_{f}}}{t_{f}^{2}} + n_{0},

Eq. (13)

z_{s_{1}}^{w} = \frac{A_{v} e^{- α v_{k} (t_{f} - T)}}{{(t_{f} - T)}^{2}} + n_{1},

where the ratio

I_{o} / v_{k}^{2}

is defined as

A_{v}

. An initialization point close enough to the actual value is required so that the Newton–Raphson algorithm converges within a limited number of iterations. This is due to the fact that the initial value of the parameter needs to be close to the true value for convergence.²² 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

Eq. (14)

γ = \frac{z_{s_{1}}^{w}}{z_{s_{0}}^{w}} = \frac{t_{f}^{2} e^{α 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, $e^{α v_{k} T}$ will be approximated by one.⁸ So with this assumption, Eq. (14) becomes

Eq. (15)

t_{f} = \frac{\sqrt{γ} T}{\sqrt{γ} - 1} .

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

Eq. (16)

A_{v} = z_{s_{0}}^{w} {(\frac{\sqrt{γ} T}{\sqrt{γ} - 1})}^{2} .

Equations (15) and (16) are then used to initialize the estimator for TTG, thus providing us with near optimal results. Next, the CRLB on the TTG estimation error is calculated. The CRLB is a lower bound on the variance of any unbiased estimate. The CRLB is generally used to determine performance of an estimator as it is referred to as the best achievable performance. The objectives of deriving the CRLB are threefold. The first one is to check if the best achievable TTG estimate is useful in enhancing the system observability, i.e., improving the target state estimates. It is shown in later sections that the TTG estimate does improve the overall performance when used as an augmented measurement to the angle information of the two IR sensors. The second one is to measure the performance of the TTG estimator and the last one is to use the lower bound as the variance of the TTG information used as an augmented measurement of the EKF for target state estimation for target ranging.

Suppose an unbiased estimate of the parameter set $θ$ 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 [{(\hat{θ} - θ)}^{2}] \geq I^{- 1} (θ)$ . The information matrix can be computed by taking the expected value of the Hessian of the log-likelihood function as follows:

Eq. (17)

I (θ) = E {{[\nabla_{θ} \log p (z_{s_{k}}^{w} | θ)]}^{T} [\nabla_{θ} \log p (z_{s_{k}}^{w} | θ)]} .

The likelihood function that will be used to calculate the information matrix can be written as

Eq. (18)

p (z_{s_{k}}^{w} | θ) = N (z_{s_{k}}^{w}; h_{k}^{w}, σ_{w}^{2}),

and with measurement vector

Z_{s}^{w}

, the likelihood function becomes

Eq. (19)

p (Z_{s}^{w} | θ) = \prod_{k = 0}^{n - 1} N (z_{s_{k}}^{w}; h_{k}^{w}, σ_{w}^{2}), = \frac{1}{{(2 π)}^{n / 2} σ_{w}^{n}} e^{- \frac{1}{2 σ_{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

Eq. (20)

h_{k}^{w} = \frac{I_{o}}{v^{2}} \frac{1}{{(t_{f} - k T)}^{2}} = \frac{A_{v}}{{(t_{f} - k T)}^{2}},

where

n

is the sample size, i.e.,

t_{f} = n T

and

k = 0, 1, 2, \dots, n - 1

. The gradient can be written as

\nabla_{θ} = [\frac{\partial}{\partial θ_{1}}, \frac{\partial}{\partial θ_{2}}], θ = {[t_{f}, A_{v}]}^{T},

where the parameter vector

θ

consists of two elements

t_{f}

and

A_{v}

. Taking the natural log of the likelihood function Eq. (19) results in

Eq. (21)

P = \log p (Z_{s}^{w} | θ), = - \log ({(2 π)}^{n / 2} σ_{w}^{n}) - \frac{1}{2 σ_{w}^{2}} \sum_{k = 0}^{n - 1} {(z_{s_{k}}^{w} - h_{k}^{w})}^{2},

and differentiating Eq. (21) with respect to

t_{f}

yields

Eq. (22)

\frac{\partial P}{\partial t_{f}} = \frac{1}{2 σ_{w}^{2}} \sum_{k = 0}^{n - 1} 2 (z_{s_{k}}^{w} - h_{k}^{w}) (\frac{\partial h_{k}^{w}}{\partial t_{f}}),

after differentiating the term

h_{k}^{w}

, given by Eq. (20), then Eq. (22) becomes

Eq. (23)

\frac{\partial P}{\partial t_{f}} = \frac{- 2 A_{v}}{σ_{w}^{2}} \sum_{k = 0}^{n - 1} \frac{(z_{s_{k}}^{w} - h_{k}^{w})}{{(t_{f} - k T)}^{3}},

similarly differentiating Eq. (21) with respect to

A_{v}

yields

Eq. (24)

\frac{\partial P}{\partial A_{v}} = \frac{1}{σ_{w}^{2}} \sum_{k = 0}^{n - 1} \frac{(z_{s_{k}}^{w} - h_{k}^{w})}{{(t_{f} - k T)}^{2}} .

By combining Eqs. (23) and (24), the gradient of the log-likelihood function with respect to the vector $θ$ can now be written as

Eq. (25)

\nabla_{θ} \log p (Z_{s}^{w} | θ) = {[\begin{matrix} \frac{- 2 A_{v}}{σ_{w}^{2}} \sum_{k = 0}^{n - 1} \frac{(z_{s_{k}}^{w} - h_{k}^{w})}{{(t_{f} - k T)}^{3}} \\ \frac{1}{σ_{w}^{2}} \sum_{k = 0}^{n - 1} \frac{(z_{s_{k}}^{w} - h_{k}^{w})}{{(t_{f} - k T)}^{2}} \end{matrix}]}^{'} .

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

Eq. (26)

I (θ) = (\begin{matrix} E {D_{11}} & E {D_{12}} \\ E {D_{21}} & E {D_{22}} \end{matrix}),

where each element of this matrix is obtained after manipulating and taking the expected value of Eq. (17). The first diagonal entry related to the TTG can be expressed as

Eq. (27)

E {D_{11}} = \frac{4 A_{v}^{2}}{σ_{w}^{2}} \sum_{k = 0}^{n - 1} \frac{1}{{(t_{f} - k T)}^{6}},

while the two nondiagonal entries take the same value and can be written as

Eq. (28)

E {D_{12}} = E {D_{21}} = \frac{- 2 A_{v}}{σ_{w}^{2}} \sum_{k = 0}^{n - 1} \frac{1}{{(t_{f} - k T)}^{5}},

and the second diagonal entry related to

A_{v}

can be written as

Eq. (29)

E {D_{22}} = \frac{1}{σ_{w}^{2}} \sum_{k = 0}^{n - 1} \frac{1}{{(t_{f} - k T)}^{4}} .

Finally, the CRLB is calculated using Eqs. (27)–(29) in Eq. (26) and then inverting the Fisher information matrix. Then

Eq. (30)

CRLB = {([P_{0 | 0}^{- 1} + I (θ)])}^{- 1},

where the first and second diagonal entries of Eq. (30) represent the CRLB of the desired parameters and

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.¹⁸ The target’s trajectory state propagates in the discrete domain by

Eq. (31)

X_{k} = F_{k} X_{k - 1} + v_{k},

where

X_{k}

is the target state vector and

v_{k}

is zero-mean white Gaussian noise with known covariance matrix

Q

. The target state vector is defined by

Eq. (32)

X = {[\begin{matrix} x & y & z & \dot{x} & \dot{y} & \dot{z} \end{matrix}]}^{'} .

The state transition matrix $F$ in Eq. (31) is defined as

Eq. (33)

F = (\begin{matrix} I_{3} & T I_{3} \\ O_{3} & I_{3} \end{matrix}), I_{3} = (\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix}),

where

O_{3}

is a

3 \times 3

zero matrix. The covariance matrix

Q

defining process noise

v_{k}

can be expressed as

Eq. (34)

Q = q . (\begin{matrix} \frac{T^{4}}{4} I_{3} & \frac{T^{3}}{2} I_{3} \\ \frac{T^{3}}{2} I_{3} & T^{2} I_{3} \end{matrix}),

where

q

is the process noise variance and

I_{3}

is the same as defined in Eq. (31). The measurement model is completely defined by Eqs. (1)–(4). The EKF update and prediction equations as given in the literature²^,¹⁸ can be expressed as

Eq. (35)

{\hat{X}}_{k / k - 1} = {F \hat{X}}_{k - 1 / k - 1},

Eq. (36)

{\hat{P}}_{k / k - 1} = {F \hat{P}}_{k - 1 / k - 1} F^{T} + Q,

whereas the estimation part can be written as

Eq. (37)

S_{k} = H_{k}^{θ} {\hat{P}}_{k / k - 1} {(H_{k}^{θ})}^{T} + R^{θ},

Eq. (38)

K_{k} = {\hat{P}}_{k / k - 1} {(H_{k}^{θ})}^{T} S_{k}^{- 1},

Eq. (39)

{\hat{X}}_{k / k} = {\hat{X}}_{k / k - 1} + K_{k} [z_{{(1,2)}_{k}}^{θ} - h_{k}^{θ} ({\hat{X}}_{k / k - 1})],

Eq. (40)

{\hat{P}}_{k / k} = (I - K_{k} H_{k}^{θ}) {\hat{P}}_{k / k - 1},

where the nonlinear function

h_{(1,2)}^{θ} ({\hat{X}}_{k / k - 1})

is a

4 \times 1

vector as defined by Eq. (3) and

H_{k}^{θ}

is the Jacobian matrix defined as

Eq. (41)

H_{k}^{θ} = (\begin{matrix} \frac{\partial η_{1_{k}}}{\partial x} & \frac{\partial η_{1_{k}}}{\partial y} & 0 & 0 & 0 & 0 \\ \frac{\partial ε_{1_{k}}}{\partial x} & \frac{\partial ε_{1_{k}}}{\partial y} & \frac{\partial ε_{1_{k}}}{\partial z} & 0 & 0 & 0 \\ \frac{\partial η_{2_{k}}}{\partial x} & \frac{\partial η_{2_{k}}}{\partial y} & 0 & 0 & 0 & 0 \\ \frac{\partial ε_{2_{k}}}{\partial x} & \frac{\partial ε_{2_{k}}}{\partial y} & \frac{\partial ε_{2_{k}}}{\partial z} & 0 & 0 & 0 \end{matrix}) .

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

\frac{\partial η_{s}}{\partial x} = \frac{y^{s} - y^{t}}{{(x^{t} - x^{s})}^{2} + {(y^{t} - y^{s})}^{2}},

\frac{\partial η_{s}}{\partial y} = \frac{x^{t} - x^{s}}{{(x^{t} - x^{s})}^{2} + {(y^{t} - y^{s})}^{2}},

\frac{\partial ε_{s}}{\partial x} = \frac{- (x^{t} - x^{s}) (z^{t} - z^{s})}{D^{2} \sqrt{{(x^{t} - x^{s})}^{2} + {(y^{t} - y^{s})}^{2}}},

\frac{\partial ε_{s}}{\partial y} = \frac{- (y^{t} - y^{s}) (z^{t} - z^{s})}{D^{2} \sqrt{{(x^{t} - x^{s})}^{2} + {(y^{t} - y^{s})}^{2}}},

\frac{\partial ε_{s}}{\partial z} = \frac{\sqrt{{(x^{t} - x^{s})}^{2} + {(y^{t} - y^{s})}^{2}}}{D^{2}},

where

D

is defined as

D = \sqrt{[{(x^{t} - x^{s})}^{2} + {(y^{t} - y^{s})}^{2} + {(z^{t} - z^{s})}^{2}]} .

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 $R^{θ}$ 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 m / 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.

Fig. 3

Target state estimation with the extended measurement vector: block diagram of the complete algorithm.

4.1.

Time-To-Go Estimation

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

Eq. (42)

\frac{\partial h_{k}^{w}}{\partial θ_{N}} = (\frac{- 2 A_{v}}{{(t_{f} - k T)}^{3}}, \frac{1}{{(t_{f} - k T)}^{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 $z_{s_{k}}^{w}$

if range less than bound then

Gather minimum required input samples.

if required samples criteria are met

calculate

θ

for start point using (15) and (16)

for number of iterations do

Calculate:

M = \sum_{k = 0}^{n - 1} {(\partial h_{k}^{w} / \partial θ_{N})}^{T} σ_{w}^{- 2} (\partial h_{k}^{w} / \partial θ_{N})

N = \sum_{k = 0}^{n - 1} {(\frac{\partial h_{k}^{w}}{\partial θ_{N}})}^{T} σ_{w}^{- 2} (z_{s_{k}}^{w} - h_{k}^{w})

Update:

{\hat{θ}}_{n}^{l + 1} = {\hat{θ}}_{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 $Z_{s_{k}}$ and hence an extended measurement matrix $H_{k}$ . The new measurement vector, which consists of five elements, i.e., target azimuth and elevation from two sensors along with TTG $(τ_{k})$ , can be expressed as

Eq. (43)

Z_{s_{k}} = {[\begin{matrix} η_{1_{k}} & ε_{1_{k}} & η_{2_{k}} & ε_{2_{k}} \end{matrix}]}^{T} .

The measurement matrix will have an additional row of partial derivatives of $(τ_{k})$ with respect to the state vector, which can be written as

Eq. (44)

H_{k}^{τ} = [\begin{matrix} \frac{\partial τ_{k}}{\partial x} & \frac{\partial τ_{k}}{\partial y} & \frac{\partial τ_{k}}{\partial z} & \frac{\partial τ_{k}}{\partial \dot{x}} & \frac{\partial τ_{k}}{\partial \dot{y}} & \frac{\partial τ_{k}}{\partial \dot{z}} \end{matrix}],

where

Eq. (45)

τ = - \frac{R}{\dot{R}} = - \frac{x^{2} + y^{2} + z^{2}}{x \dot{x} + y \dot{y} + z \dot{z}},

with

R

being the range and

\dot{R}

the range rate. These are calculated based on the target state prediction and not based on the actual position. The new measurement matrix is obtained by combining Eq. (44) with Eq. (41) and can be written as

Eq. (46)

H_{k} = [\begin{matrix} H_{k}^{θ} \\ H_{k}^{τ} \end{matrix}] .

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 m / 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 m$ . As illustrated in Fig. 4, the variance of the range with a baseline length of $b = 1 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 m$ .

Fig. 4

Variance, $σ_{R_{k}}^{2}$ , of range measurements with different baseline lengths.

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,²³^,²⁴ 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 $σ_{η} = 5.7 mrad$ and $σ_{ε} = 5.7 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 s$ are given by Figs. 5 and 6.

Fig. 5

Square root of Cramér–Rao lower bound (CRLB) for $t_{f}$ .

Fig. 6

Square root CRLB for parameter $A_{v}$ .

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, $R^{R}$ , and are quantified by standard deviations, $σ_{r} = 1 m$ , $σ_{ε} = 26.4 mrad$ , and $σ_{η} = 26.4 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 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.

Fig. 7

(a) Root mean square error (RMSE) in time-to-go (TTG) estimate as a function of target range (500 to 100 m) for $SNR = 25 dB$ . (b) RMSEs in TTG estimates and the CRLBs for various SNRs at target range of 100 m.

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.

Fig. 8

RMSE in position estimate.

Fig. 9

RMSE in velocity estimate.

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 $TTG = 0.4 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, $σ_{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 $TTG = 0.25 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

1.

J. H. Wehling, “Multifunction millimeter-wave systems for armored vehicle application,” IEEE Trans. Microwave Theory Tech., 53 1021 –1025 (2005). http://dx.doi.org/10.1109/TMTT.2005.843504 IETMAB 0018-9480 Google Scholar

2.

S. Blackman and R. Popoli, Design and Analysis of Modern Tracking Systems, Artech House, Boston, Massachusetts (1999). Google Scholar

3.

P. Wu et al., “Tracking algorithm with radar and infrared sensors using a novel adaptive grid interacting multiple model,” IET Sci. Meas. Technol., 8 270 –276 (2014). http://dx.doi.org/10.1049/iet-smt.2013.0020 ISMTCT 1751-8822 Google Scholar

4.

R. Mobus and U. Kolbe, “Multi-target multi-object tracking, sensor fusion of radar and infrared,” in IEEE Intelligent Vehicles Symposium, 732 –737 (2004). http://dx.doi.org/10.1109/IVS.2004.1336475 Google Scholar

5.

M. de Visser et al., “Passive ranging using an infrared search and track sensor,” Opt. Eng., 45 (2), 026402 (2006). http://dx.doi.org/10.1117/1.2173948 OPEGAR 0091-3286 Google Scholar

6.

F. Dufour and M. Mariton, “Target tracking with two passive infrared non-imaging sensors,” IEEE Trans. Aerosp. Electron. Syst., 27 (4), 725 –739 (1991). http://dx.doi.org/10.1109/7.85047 IEARAX 0018-9251 Google Scholar

7.

Z. Djurovic, B. Kovacevic and G. Dikic, “Tracking algorithm with radar and infrared sensors using a novel adaptive grid interacting multiple model,” IET Signal Process., 3 (3), 177 –188 (2009). http://dx.doi.org/10.1049/iet-spr.2008.0068 1751-9675 Google Scholar

8.

S. Kim, S. G. Sun and K. T. Kim, “Analysis of infrared signature variation and robust filter-based supersonic target detection,” Sci. World J., 2014 (3), 17 (2014). http://dx.doi.org/10.1155/2014/140930 1597-6343 Google Scholar

9.

M. Ito, S. Tsujimichi and Y. Kosuge, “Tracking a three-dimensional moving target with distributed passive sensors using extended Kalman filter,” Electron. Commun. Jpn., 84 (7), 74 –85 (2001). http://dx.doi.org/10.1002/ecja.1026 ECJCED 8756-6621 Google Scholar

10.

B. Ristic and M. S. Arulampalam, “Tracking a manoeuvring target using angle-only measurements: algorithms and performance,” Signal Process., 83 (6), 1223 –1238 (2003). http://dx.doi.org/10.1016/S0165-1684(03)00042-2 SPRODR 0165-1684 Google Scholar

11.

P. Wu and X. Li, “Passive multi-sensor maneuvering target tracking based on UKF-IMM algorithm,” in WASE Int. Conf. on Information Engineering (ICIE’09), 135 –38 (2009). http://dx.doi.org/10.1109/ICIE.2009.182 Google Scholar

12.

J. Zhang et al., “Infrared target tracking using multisensor data fusion,” Proc. SPIE, 6835 68351W (2008). http://dx.doi.org/10.1117/12.757379 PSISDG 0277-786X Google Scholar

13.

T. A. Camus, “Calculating time-to-collision with real-time optical flow,” Proc. SPIE, 2308 661 –670 (1994). http://dx.doi.org/10.1117/12.186009 PSISDG 0277-786X Google Scholar

14.

N. Ancona and T. Poggio, “Optical flow from 1-D correlation: application to a simple time-to-crash detector,” Int. J. Comput. Vision, 14 (2), 131 –146 (1995). http://dx.doi.org/10.1007/BF01418979 IJCVEQ 0920-5691 Google Scholar

15.

G. Alenya, A. Nègre and J. L. Crowley, “A comparison of three methods for measure of time to contact,” 4565 –4570 IEEE(2009). http://dx.doi.org/10.1109/IROS.2009.5354024 Google Scholar

16.

R. Cipolla and A. Blake, “Surface orientation and time to contact from image divergence and deformation,” in Computer Vision-ECCV, 187 –202 (1992). http://dx.doi.org/10.1007/3-540-55426-2-21 Google Scholar

17.

A. Fernández-Caballero et al., “Optical flow or image subtraction in human detection from infrared camera on mobile robot,” Rob. Auton. Syst., 58 (12), 1273 –1281 (2010). http://dx.doi.org/10.1016/j.robot.2010.06.002 RASOEJ 0921-8890 Google Scholar

18.

B. Ristic, S. Arulampalam and N. Gordon, Beyond the Kalman Filter: Particle Filters for Tracking Applications, 685 Artech House, Boston (2004). Google Scholar

19.

G. Waldman and J. Wootton, Electro-Optical Systems Performance Modeling, Artech House, Boston (1993). Google Scholar

20.

Y. Bar-Shalom and X. R. Li, Multitarget–Multisensor Tracking: Principles and Techniques, YBS PublishingUrbana, Illinois,1995). Google Scholar

21.

A. Cherchar, A. Belouchrani and T. Chonavel, “A new approach to maneuvring target tracking in passive multisensor environment,” (2010). Google Scholar

22.

B. Bradie, A Friendly Introduction to Numerical Analysis, Pearson Prentice Hall, Upper Saddle River, New Jersey (2006). Google Scholar

23.

THALES, “Compact high performance megapixel thermal imager,” (2014) https://www.thalesgroup.com/sites/default/files/asset/document/1.1 - Cath MP Datasheet.pdf October ). 2014). Google Scholar

24.

FLIR Systems, “FLIR HDC, long-range thermal imaging cameras,” (2014) http://www.flirmedia.com/MMC/CVS/Comm_sec/SS_0030_EN.pdf October ). 2014). Google Scholar

25.

D. Lerro and Y. Bar-Shalom, “Tracking with debiased consistent converted measurements versus EKF,” IEEE Trans. Aerosp. Electron. Syst., 29 (3), 1015 –1022 (1993). http://dx.doi.org/10.1109/7.220948 IEARAX 0018-9251 Google Scholar

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.

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.

Citation Download Citation

Ihsan Ullah, Taek Lyul Song, and Thia Kirubarajan "Active vehicle protection using angle and time-to-go information from high-resolution infrared sensors," Optical Engineering 54(5), 053110 (13 May 2015). https://doi.org/10.1117/1.OE.54.5.053110

Published: 13 May 2015

Access the abstract

JOURNAL ARTICLE
10 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY