Bearings-Only Target Motion Analysis with Acoustic Propagation Models of Uncertain Fidelity

The main theoretical contribution of the work presented here is the insight that partial state observability problems can be treated as “missing data” problems for which the expectationmaximization (EM) method is a potentially useful technique for deriving new state estimation algorithms. The missing data model can be applied to general partial state observability problems; however, attention is focused here exclusively on its application to the passive sonar target motion analysis (TMA) problem. The TMA application is somewhat narrowly specialized, but it exposes the main ideas and issues encountered in applications of missing data models to partial state observability. A bearing measurement does not fully characterize target location; that is, in the parlance of EM, the corresponding range measurement is “missing.” The missing data in the bearings-only TMA problem are the missing range measurements. Modeling the joint probability density function (pdf) of the observed measurements and their missing data in a manner suitable for application of the EM method is the crux of the missing data model. The joint density is chosen in an application specific manner, and it ultimately determines the utility of the state estimation algorithm derived via EM. For the joint pdf developed here for the bearings-only TMA problem, application of the EM method yields an iteratively reweighted linear least-squares algorithm. The notation needed to derive this algorithm is burdensome, but the algorithm itself is simple, intuitive, and strikingly different from the traditional nonlinear least-squares algorithm derived via the maximum likelihood (ML) method. The classical bearings-only TMA problem is a special case of a more general TMA problem that is called herein the augmented bearings-only TMA problem. The augmented bearings-only TMA problem is defined here as the problem that arises when bearing measurements are augmented with measurements of the received signal-to-noise ratio (SNR). A combined acoustic propagation and sensor (CAPS) performance prediction model is assumed given. The fidelity of available ocean acoustic propagation models is often not completely satisfactory; that is, accurate models of the environment may yield predictions that are mismatched to the real world because of the lack of adequate or current environmental data. Nonetheless, it may be possible to exploit such models for TMA purposes. One contribution of our work is a novel mathematical treatment, first reported in [1, 2], of the environmental model fidelity issue. There are two reasons to generalize the classical bearings-only TMA problem. One reason is that the augmented problem is potentially useful in


I. INTRODUCTION
The main theoretical contribution of the work presented here is the insight that partial state observability problems can be treated as "missing data" problems for which the expectationmaximization (EM) method is a potentially useful technique for deriving new state estimation algorithms. The missing data model can be applied to general partial state observability problems; however, attention is focused here exclusively on its application to the passive sonar target motion analysis (TMA) problem. The TMA application is somewhat narrowly specialized, but it exposes the main ideas and issues encountered in applications of missing data models to partial state observability.
A bearing measurement does not fully characterize target location; that is, in the parlance of EM, the corresponding range measurement is "missing." The missing data in the bearings-only TMA problem are the missing range measurements. Modeling the joint probability density function (pdf) of the observed measurements and their missing data in a manner suitable for application of the EM method is the crux of the missing data model. The joint density is chosen in an application specific manner, and it ultimately determines the utility of the state estimation algorithm derived via EM. For the joint pdf developed here for the bearings-only TMA problem, application of the EM method yields an iteratively reweighted linear least-squares algorithm. The notation needed to derive this algorithm is burdensome, but the algorithm itself is simple, intuitive, and strikingly different from the traditional nonlinear least-squares algorithm derived via the maximum likelihood (ML) method.
The classical bearings-only TMA problem is a special case of a more general TMA problem that is called herein the augmented bearings-only TMA problem. The augmented bearings-only TMA problem is defined here as the problem that arises when bearing measurements are augmented with measurements of the received signal-to-noise ratio (SNR). A combined acoustic propagation and sensor (CAPS) performance prediction model is assumed given. The fidelity of available ocean acoustic propagation models is often not completely satisfactory; that is, accurate models of the environment may yield predictions that are mismatched to the real world because of the lack of adequate or current environmental data. Nonetheless, it may be possible to exploit such models for TMA purposes. One contribution of our work is a novel mathematical treatment, first reported in [1,2], of the environmental model fidelity issue.
There are two reasons to generalize the classical bearings-only TMA problem. One reason is that the augmented problem is potentially useful in applications. The other reason is purely technical; that is, the missing data model leads to an integral representation involving a parametric target density called the geometric kernel (see Section IIB below). The function multiplying the geometric kernel is identically equal to one for classical bearings-only TMA and, consequently, important features of the missing data model are not fully exposed. The augmented TMA problem, however, leads to nontrivial multipliers, and the resulting integral representation can be interpreted as an integral transform of the CAPS density.
The TMA formulation presented here is statistically unconventional in that it assumes the measurements condition a sequence of independent, predetection "empirical" target random variables. Multiplying the empirical target pdfs gives the joint density on the empirical target states. Empirical maximum a posteriori (EMAP) target state estimates are obtained by maximizing this joint pdf for a specified class of parametric target motion models. EMAP estimates are computed-without taking the gradient of the CAPS model-using an iteratively reweighted, linear least-squares algorithm derived via the EM method. Consequently this algorithm also solves the classical (nonaugmented) bearings-only TMA problem [3]. In contrast to the ML approach, which generally requires computing gradients of the CAPS model, the EMAP approach requires computing integrals of the CAPS model with respect to the geometric kernel. Adjusting the down-range variance of the geometric kernel is an intuitive heuristic method for compensating for CAPS model mismatch. The methods proposed here are directly applicable to other passive sonar problems, e.g. TMA using cone angle measurements from linear arrays [4]; however these problems are outside the scope of this work.
The EMAP formulation of the augmented bearings-only TMA problem is presented in Section II. The CAPS model pdf is discussed in Section III, and the EMAP and ML approaches are compared in Section IV. EMAP estimation is discussed in Section V, and a constant velocity target example is presented in Section VI. Concluding remarks are given in Section VII. Appendix A shows that the asymptotic limit of the EMAP formulation is equivalent to the ML formulation, and Appendix B gives the details of the derivation of the EMAP algorithm using the method of EM.

II. EMAP FORMULATION
Let X = (x, y, s) denote target location (x, y) and source (or signal) level s radiated in the azimuthal direction°, referred to as the aspect angle, measured counterclockwise from the bow of the target at an arbitrary, but fixed, time (see Fig. 1). Radiated source level may depend on elevation angle, but for simplicity, only azimuthal variation is considered here. The array signal processor estimates an azimuthal arrival angle, measured in the xy-plane counterclockwise from the positive x-axis, and an SNR, referenced to the array beamformer output. These estimates of bearing and SNR are referred to as measurements here. Let Z = (µ, ») denote the measured bearing µ and measured SNR » of the target when sensor location is ³ X = ( ³ x, ³ y). Sensor location ³ X is taken here to be a known parameter, but can be treated as a random variable if desired.

A. Empirical Target PDF
Traditional statistical models for TMA are postdetection models; that is, they assume a priori that the measurements belong to the same target. Postdetection TMA assumes that measurements are independent if they are conditioned on the target. ML target state estimators thus answer the question "Given data generated from a target track, which parameterized track best fits the data?" In contrast, the EMAP target state estimators proposed here are joint detection-estimation methods which seek to answer the alternative question "Does a target track of the specified parametric form fit the data?" In the EMAP approach, measurements are statistically independent, but parametrically tied.
Let Z´fZ n g N n=1´f µ n , » n g N n=1 denote a sequence of independent sensor measurements obtained at sensor locations ³ X´f ³ X n g N n=1´f ³ x n , ³ y n g N n=1 and at denote so-called empirical random variables associated with potential target locations and source levels. It is assumed that location (x − n , y − n ) is independent of source level s − n for all n. The random variable pairs f(Z n , X − n )g N n=1 are statistically independent because measurements are not specified a priori to belong to the same target. Independence and Bayes Theorem imply that where the joint empirical target pdf is and where is a data-dependent constant independent of X − . The joint density (2) is the basis for EMAP estimation. In the following subsection, we show that the conditional density of X − n has the form n¸¯· ³ x n + r n cos µ n ³ y n + r n sin µ n¸, r 2 n § n (µ n ) ¶ r n dr n (4) where N(¢ j ¹, C) denotes the Gaussian pdf with mean vector ¹ and covariance matrix C. The first density in the integrand of (4) is the CAPS density of the SNR measurement » n , and the second density is the pdf of empirical target location conditioned on the measured bearing µ n and the unobserved range r n . EMAP state estimates are obtained by substituting a parametric target motion model into (4), substituting this result into (2), and maximizing with respect to the target parameters (see Section VA). If the measured SNR » n is not available or deemed too unreliable to be useful, the CAPS density can be eliminated by marginalizing (4) over f» n g N n=1 , thereby obtaining a new integral formulation of the classical bearings-only TMA problem (see Section IV and [3]). In any event, the joint density (4) is a product of integrals over unobserved range measurements and, hence, amenable to treatment by the missing data methods of EM.

B. Integral Representation of Empirical Target PDF
Full target state is not observable from a single bearing measurement µ n . Consequently the dummy random variable r n is introduced to model the missing range measurement, and the conditional pdf of the empirical random variable X − n = (x − n , y − n , s − n ) is expressed as the marginal over r n of the pdf of X − n conditioned on the "complete" data (r n , µ n , » n ). Using Bayes Theorem, the marginal density is written £ p r n s − n jZ n (r n , s − n j Z n ; ³ X n )r n dr n where the Jacobian r n is needed because (r n , µ n ) is a polar coordinate system (with area differential r n dr n dµ n ). Recalling that location (x − n , y − n ) is independent of source level s − n , it follows that Several applications of Bayes Theorem gives the identity p r n s − n jZ n (r n , s − n j Z n ; ³ X n ) = p » n jr n µ n s − n (» n j r n , µ n , s − n ; ³ X n ) £ p r n jµ n (r n j µ n ; ³ X n ) p » n jµ n (» n j µ n ; ³ X n ) p s − n jr n µ n (s − n j r n , µ n ; ³ X n ): (7) Because the location measurements (r n , µ n ) do not constrain the velocity (or aspect) of the empirical target X − n , the empirical target source level s − n is independent of the measurements r n and µ n , so that p s − n jr n µ n (s − n j r n , µ n ; ³ X n )´p s − n (s − n ; ³ X n ): In the case where target source level s n is known and independent of the azimuthal angle°n, the pdf (8) is simply where S 0 is the target source level and ±(¢) denotes the Dirac delta function. Substituting (6)-(9) into (5) and marginalizing over empirical target source level s − n gives £ p x − n y − n jrnZn (x − n , y − n j r n , Z n ; ³ X n )p rnjµn (r n j µ n ; ³ X n )r n dr n (10) where ® n is a data-dependent constant.
In the general case where target source level s n is a known function of the azimuthal angle°n, the pdf (8) is written where S(°n) is the target source level (within the sensor bandwidth), and the aspect angle°n is parametrically determined by the target motion model (see Section VA, especially (22)). Target source level may also be speed dependent, but such dependence is not treated here. The densities (9) and (11) can be replaced by other pdfs to obtain fully stochastic models of source level; however, the deterministic forms (9) and (11) are used throughout this work. Each density in the integrand of (10) has a meaningful physical interpretation, as discussed in the rest of this section and in Section III, and illustrated in Section IV. 1) Geometric Kernel Density: The second density in the integrand of the empirical target pdf (10), called the geometric kernel, is a density on empirical target location. For the bearings-only TMA problem, the geometric kernel is assumed to be a bivariate Gaussian whose mean vector and covariance matrix are determined by the conditioning variables. Purely geometrical considerations (see [1,Appendix A,(A.15) and (A. 16)]) show that this Gaussian kernel must have the form n¸¯· ³ x n + r n cos µ n ³ y n + r n sin µ n¸, r 2 n § n (µ n ) ¶ : The covariance matrix § n (µ n ) in (12) may be factored in the form where U(µ) = · cos µ sin µ ¡ sin µ cos µi s a rotation matrix based on the measured bearing µ, and ¾ 2 n and · 2 n are the range-normalized cross-range and down-range variances of the kernel density, respectively.
The range-normalized cross-range variance ¾ 2 n is determined by sensor signal processing considerations (e.g., the time-bandwidth product of the received signal). It may be constant, but in general it may depend parametrically on bearing because of beamwidth equalization issues (i.e., some beams may be narrower than others), and on true (not measured) SNR. In practice, the variance ¾ 2 n can be specified as a function of the measured SNR. In the general case of aspect-angle dependent source level, the parametric form of ¾ 2 n´¾ 2 n (S(°n)) is assumed given. The range-normalized down-range variance · 2 n is a more difficult parameter to understand than ¾ 2 n because it is the variance of the missing range variable r n . As shown in Appendix A, the ML formulation of the augmented bearings-only TMA problem is obtained from the EMAP formulation in the limit as · n ! 0. In this paper · n is taken to be a free parameter. As discussed in Section III, · n is selected to compensate for CAPS model mismatch and to accelerate algorithm convergence.
2) Sensor Range pdf: The third density in the integrand of the empirical target pdf (10) specifies the sensor measurement window for the missing range measurement r n . It is derived by assuming that the a priori joint density of the measurement pair (r n , µ n ) is uniform over an annular feasible region D( ³ X n ) of the xy-plane whose inner and outer radii are given by the radial functions r (n) min´rmin (µ n , ³ X n ) and r (n) max´rmax (µ n , ³ X n ), respectively. The outer radius may be interpreted as the maximum range at which signals of specified source level are detected (with specified probability P d ). Similarly, the inner radius may be interpreted as the near-field limit of the sensor, i.e., the minimum range at which the sensor's beamformer reliably estimates bearings. In polar coordinates, the joint density of (r n , µ n ) is given by where p r n µ n (r n , µ n ; ³ X n )r n dr n dµ n = 1: Conditioning on bearing gives, using Bayes Theorem,

III. CAPS MODEL PDF
The first density in the integrand of the empirical target pdf (10) is the CAPS density of measured SNR » n conditioned on measured bearing µ n , which may in general depend on the absolute location of the sensor in the ocean. As a simple example, consider the following CAPS density for known source level, independent of time, bearing, and sensor location: (16) where´2 is the variance and a¸0 is a scalar quantity that depends on source level. A contour plot of this pdf as a function of target location (x, y) is shown in Fig. 2(a) for » = 33:08,´= 10:0, and a = 330:8 km. As is evident from Fig. 2(a), the peak of the pdf (16) occurs at r = 10:0 km in this case. It is shown in Section VI, the CAPS density (16) is consistent with propagation loss due to cylindrical spreading which, in decibels, is given by ¡10 log r (see propagation loss curve in Fig. 3(a)). Consider next the CAPS density where´2 is the variance and a, b, ¶, and ¹ are all positive scalar quantities. A contour plot of this pdf as a function of target location (x, y) is shown in   (17). The propagation loss function for this CAPS density is plotted, in decibels, in Fig. 3(b). The propagation loss curve in Fig. 3(b) is essentially identical to that in Fig. 3(a), but with a convergence zone at 25 km (see [5, sect. 6.3]).
If a deterministic CAPS model is available, it may be useful to treat it as the mean of a Gaussian pdf with suitably specified range-and bearing-dependent variances (see example in Section VI); however, such statistical models must be used with caution. An alternative approach is to use techniques such as those in [6].
The range-normalized down-range variance · 2 n of the geometric kernel provides a new and useful parameter that can be used to compensate for mismatch between the CAPS model and the real world. In effect, the CAPS model prediction is averaged over a sliding Gaussian window whose size, both down-range and cross-range, increases linearly with range. Thus, using a large value of · n decreases the CAPS model influence on TMA estimates and, conversely, using a small value increases its influence. Thus, the selection of · n depends on the degree of confidence one has in the environmental model. As · n ! 0, the kernel density approaches the product of a Dirac delta function down-range and a Gaussian density cross-range, effectively sampling the CAPS density at the kernel mean. Appendix A shows that, in the limit as the down-range length · n of the averaging window goes to zero, the EMAP approach is equivalent to the ML approach, which is the optimal approach when CAPS model predictions are assumed exact.
For sufficiently small non-zero values of · n , the EMAP algorithm will experience slow convergence rates. This convergence property is explained by results in [7,8] which show that the convergence rate of an EM algorithm is inversely proportional to the information in the missing data. Numerical experience also supports this observation. In those (highly atypical) applications in which the CAPS model is known to have such good fidelity that · n must be chosen very small, starting with a large value of · n and using a succession of progressively smaller values significantly improves convergence rates in the EMAP algorithm (see example in Section VI).

IV. COMPARISON OF EMPIRICAL TARGET PDF AND MEASUREMENT LIKELIHOOD FUNCTION
It is worthwhile to examine the empirical target pdf (4) in the special case when the measured SNR is unavailable. Integrating (4) over » n , dropping the subscript n and superscript −, and omitting the proportionality factor gives the empirical target pdf in the form A contour plot of (18) is shown in Fig. 4(a) for ³ x = ³ y = 0, r min = 500 m, r max = 30 km, µ = 0, ¾ = 0:0175 (= 1 ± ), and · = 0:0873 (= 5¾). In essence, Fig. 4(a) depicts the posterior density on empirical target location given one bearing measurement. Ignoring edge effects, inside the annulus r min < k(x ¡ ³ x, y ¡ ³ y)k < r max , the level curves are straight lines intersecting at x ¡ ³ x = y ¡ ³ y = 0. Moreover, a circular cut at a fixed range r within the annulus is a nearly Gaussian shape with standard deviation r¾. For comparison, the likelihood function for a single bearing measurement, is shown in Fig. 4(b) for the same values of ³ x, ³ y, µ, and ¾. Clearly, the empirical target pdf (18) and the likelihood function (19) are essentially identical up to a proportionality constant, that is, Asymptotic results that rigorously show the relationship between (18) and (19) as · ! 0 are given in Appendix A.
Consider next the CAPS density given by (17). A contour plot of the empirical target pdf (18) with this CAPS density inserted into the integrand is shown in Fig. 5(a) for the same values of »,´, a, b, ¶, and ¹ used in Section III. Fig. 5(a) depicts the posterior density on empirical target location given one bearing and SNR measurement. For comparison, the likelihood function for the pair (µ, ») is shown in Fig. 5(b). Again, ignoring edge effects, the likelihood function is closely approximated by the empirical target pdf. Edge effects can be eliminated entirely by using a larger value of r max or smaller values of ·.
Appendix B of [1] shows that the down-range marginal density of the empirical target pdf (4) is closely approximated by the Gaussian distribution, provided the bearing measurement standard deviation ¾ n is small, say on the order of several degrees or less. This analytical result is important because it shows (with appropriate use of diffuse priors) that the classical bearings-only TMA formulation is recovered from (4) via down-range marginalization.
Discretizing the empirical target pdf (18) gives a finite Gaussian mixture approximation of the kind used by Kronhamn in [9][10][11] in conjunction with a multihypothesis Kalman filtering approach. Although Kronhamn's methods differ significantly from those used here, his insight that a Gaussian sum can approximate the likelihood function (19) is strikingly similar to the continuous mixture representation (18).

A. Algorithm Formulation
Let X´fX n g N n=1´f x n , y n , s n g N n=1 denote the sequence of target states at the measurement times ft n g N n=1 . A deterministic target motion model is assumed specified, so that X n is a function of t n and the target parameterization. The standard target motion model for bearings-only TMA is constant velocity, so for the end-point parameterization¸= [x 1 , y 1 , x N , y N ] T , · x n y n¸= H(t n )¸ (20) where the time-dependent measurement matrix H(t) is given by If there are no ocean currents, the bow of the target points in the direction of its motion, and simple geometrical considerations in the azimuthal plane give the source level reference angle°n (see Fig. 1). Substituting the velocity implicit ini nto the aspect-angle dependent source level function S(°n)´S(°n(¸)) gives s n = S(°n(¸)). Target state is therefore fully parameterized by¸, so that X = X(¸) = fx(¸), y(¸), s(¸)g.
To simplify notation, let ¹ n (r n , µ n )´· ³ x n + r n cos µ n ³ y n + r n sin µ n¸ (  23) where the over-set arrow is used to emphasize that ¹ n is the target location vector (with respect to the origin) given by the bearing and missing-range measurements (µ n , r n ). Substituting the target model parameterization (20)-(23) into the empirical target pdf (4), and then substituting the result into the joint pdf (2) marginalized over empirical source level gives the joint empirical target pdf as a product of integrals, i.e., min p » n jr n µ n s − n (» n j r n , µ n , S(°n(¸)); ³ X n ) £ N(H(t n )¸j1 n (r n , µ n ), r 2 n § n (µ n ))r n dr n , (24) where the range-normalized cross-range and down-range variances of the geometric kernel are ¾ 2 n´¾ 2 (S(°n(¸))) and · 2 n . The dependence of ¾ 2 n on the target parameter¸arises via target source level dependence on aspect angle.
The joint empirical target pdf (24) yields the EMAP target parameter estimatȩ The standard necessary conditions for the EMAP estimate (25) are found by setting the gradient of (24) with respect to the target parameters¸to zero. It is a remarkable fact that these conditions do not require computing the gradient of the CAPS model if target source level is independent of aspect angle; that is, if S(°) is independent of°.

B. Algorithm Statement
The EMAP target parameter estimation algorithm is derived in Appendix B using the EM method. The algorithm is derived for constant velocity target motion; however, the derivation is easily extended to general models in which the target motion parameters appear linearly. This section presents the main results from Appendix B, and lists the steps of the EMAP algorithm.
The objective of the E-step is to define the so-called auxiliary function ª (¸) of the EM method and to simplify it if possible. The auxiliary function depends on two sets of target parameters,¸0 = [x 0 1 , y 0 1 , x 0 N , y 0 N ] T and¸= [x 1 , y 1 , x N , y N ] T , where¸0 is an initial estimate and¸is arbitrary. Defining the quadratic form Q(r n , µ n ;¸) = [H(t n )¸¡1 n (r n , µ n )] T £ § ¡1 n (µ n )[H(t n )¸¡1 n (r n , µ n )] (26) the terms of the auxiliary function that are functions of¸are min p miss(n) (r n j Z n ; ³ X n ,¸0)Q(r n , µ n ;¸) dr n r n (28) is a weighted mean-squared error, min p miss(n) (r n j Z n ; ³ X n ,¸0) £ log[p »njrnµns − n (» n j r n , µ n , S(°n(¸)); ³ X n )]r n dr n (29) is a CAPS model contribution due to target source level dependence on aspect angle, and min p miss(n) (r n j Z n ; ³ X n ,¸0) £ log[¾(S(°n(¸)))]r n dr n is a variance contribution due to target source level dependence on aspect angle. The missing-range conditional pdfs in (28)-(30) are given by p miss(n) (r n j Z n ; ³ X n ,¸) = Ã miss(n) (r n j Z n ; ³ X n ,¸) R r (n) max r (n) min Ã miss(n) (½ n j Z n ; ³ X n ,¸)½ n d½ n where Ã miss(n) (r n j Z n ; ³ X n ,¸) = p » n jr n µ n s − n (» n j r n , µ n , S(°n(¸)); ³ X n ) £ N(H(t n )¸j1 n (r n , µ n ), r 2 n § n (µ n )) (32) are the unnormalized missing-range pdfs. The missing-range densities (31) are nonnegative, so ª MSE (¸) is a nonnegative quadratic function of¸. If target source level is independent of aspect angle, the term ª CAPS (¸) is omitted from the auxiliary function (27) because it is independent of¸. Similarly, the term ª VAR (¸) is omitted from (27) if the cross-range variance of the geometric kernel is independent of the target motion model. The objective of the M-step is to maximize the auxiliary function as a function of the target parameters¸. Consider first the important special case where the auxiliary function comprises only the weighted MSE term ª MSE (¸). The missing-range densities in (28) depend on the initial parameter vectoŗ 0 , but not on¸; hence, setting the gradient of ª (¸) with respect to¸to zero and solving for¸gives the updated parameter estimatȩ +´a rg max ª (¸)´arg miņ ª MSE (¸) where A n (Z n ) is the 4 £ 4 matrix given by min p miss(n) (r n j Z n ; ³ X n ,¸0) dr n r n (34) and b n (Z n ) is the vector of length 4 given by min p miss(n) (r n j Z n ; ³ X n ,¸0)1 n (r n , µ n ) dr n r n : The matrix in (34) has rank at most 2, so the matrix whose inverse is required in (33) attains full rank if, and only if, the target is observable (in the statistical sense) from the measurements. Observability questions are widely discussed in bearings-only TMA problems (see, e.g., [12]), but lie outside the intended scope of this work. The EMAP algorithm is an iteratively reweighted, linear least-squares algorithm for the important special case when only the weighted MSE term ª MSE (¸) is retained. Let S 0 denote the given target source level, independent of aspect angle. Let ¾ n and · n for n = 1,:::, N be specified. Explicitly, the algorithm for this special case takes the following recursive form.

1) Initialize the target parameter vectoŗ
N ] T , and set k = 0.

3) Using the integrals (37), compute the 4 £ 4 matrix
and the length-4 vector In practice, linear least-squares problems such as (33) are best solved by robust methods such as QR factorization or singular value decomposition, rather than explicit solution of the normal equations. Also, in general, the integrals (37) are not integrable in closed form; however, they are readily evaluated by well known and straightforward numerical techniques.
Convergence tests for the EMAP algorithm are based on the rate of increase of the joint empirical target pdf (24), which is guaranteed by the method of EM to increase monotonically at each iteration. From (24), using the integral forms (37) and omitting irrelevant multiplicative constants, the joint empirical target pdf at the kth iteration is given by Alternatively, it might be useful in practice to base a convergence test on changes in the missing range pdfs (31) (see example in Section VI).

C. Other Considerations
The EMAP algorithm is easily generalized for unknown source level by estimating it between EM iterations. This extension of the EM method is referred to as the generalized EM (GEM) method, and is discussed in [13]. Source level is estimated in the GEM framework by conducting a one-dimensional search on source level to increase the value of the auxiliary function between EM iterations. Though this approach presents no theoretical difficulties, the EMAP algorithm is more complex than a simple iteratively reweighted, linear least-squares algorithm in this case; source level estimation is important in practice but is not pursued further here.
The M-step cannot neglect the terms ª CAPS (¸) and ª VAR (¸) in the general case when target source level depends on aspect angle. Solving for the updated parameter vector in the M-step requires solving the following problem: (41) The optimization problem (41) is a nonlinearly penalized linear least-squares problem and must be solved by iterative methods such as Gauss-Newton. Unfortunately, evaluating the gradient of the objective function in (41) requires taking the derivative of the CAPS model. The EMAP algorithm for the general case, therefore, requires solving a difficult numerical problem at each iteration, and it is thus comparable in computational complexity to algorithms based on the standard ML approach. Whether the penalized least-squares formulation (41) can be exploited to advantage in the general case is an open question.

A. CAPS Model Density
It is assumed in the following example that target source level is a known constant, independent of aspect angle, and that the target signal is of sufficiently low frequency so that propagation loss due to absorption is negligible. A homogeneous, isovelocity (1500 m/s) ocean is assumed with direct-path propagation. An isotropic ambient noise field is also assumed.
The stochastic model used in this example for the SNR measurement »´»(r, t) is where»(r) is the mean SNR derived from a deterministic CAPS model and q(t) is a stationary, zero-mean, white-Gaussian-noise process with variance´2. The CAPS model mean»(r) for this simple example is given by the passive sonar equation, which, in decibels, relates signal-to-noise ratio SNR to target source level SL, propagation loss PL, ambient noise level NL, and array directivity index DI (see [14, sect. 15.8] for details). In decibels, the passive sonar equation is written as where the combination SL ¡ PL + DI is the received signal level. For this example, the notional target is assumed to have a source spectrum level SSL of 100 dB//1 ¹Pa near 1 kHz (see [5, ch. 10] for nominal values), so that SL´SSL + 10 log ¢f (44) for broadband detection, where ¢f is the receiver bandwidth in Hertz. Propagation loss is assumed due to cylindrical spreading. Thus, for a point source, the acoustic intensity (which is proportional to the square of the pressure) falls off as r, so that in decibels The propagation loss curve given by (45) is plotted in Fig. 3(a). Noise level in terms of the noise spectrum level NSL is given by NL´NSL + 10 log ¢f:   Table 3.2]). The cone angles estimated by the line array are interpreted as azimuthal bearings, with the additional assumptions that the target and the sensor are at equal depths, and that the port/starboard ambiguity associated with each cone angle is resolved. In this example, r is the distance between the acoustic center of the array and the target.
Substituting (44)-(46) into (43) gives the mean SNR, in decibels, as 10 log»(r) = SSL ¡ 10 log r + DI ¡ NSL: Substituting (47) into (42) and writing the result as a pdf gives p » n jr n s − n (» n j r n , SL)´N(» n j»(r n ),´2) for the CAPS density of the measured SNR » n . The pdf (48) is plotted in Fig. 2(a) as a function of target location (x n , y n ) for the values of SSL, DI, NSL, andú sed in this example. Fig. 6 shows a two-leg constant velocity target scenario in which the sensor moves on a fixed course of 5 ± at a speed of 5 m/s on the first leg, followed by a fixed course of 95 ± at the same speed on the second leg after a 2 min simulated maneuver (recall that all angles are measured counterclockwise from the positive x-axis). The target is initially 10 km due east of the sensor, and moves with a constant velocity of 2.5 m/s on westerly course. The sensor generates bearing and SNR measurements with standard deviations ¾ = 1 ± and´= 10, respectively, at 1 min intervals on both 10 min legs.

B. Constant Velocity Target Scenario
The EMAP algorithm given in Section V, with the CAPS density defined by (47) and (48), is used to obtain EMAP estimates of the target parameters¸. Gauss-Legendre quadrature with 50 sample ranges is used to calculate the integrals (37), where the minimum and maximum ranges are taken to be 500 m and 30 km, respectively. A dimensionless down-range standard deviation at r = 1 m of · = 0:0873 is used; the dimensionless cross-range (bearing) standard deviation at r = 1 m corresponding to ¾ = 1 ± is ¾ = 0:0175. These values for · and ¾ translate roughly to down-range and cross-range standard deviations for a target r = 10 km away of · r´· r = 873 m and ¾ r´¾ r = 175 m, respectively. ML estimates of¸are obtained by maximizing the likelihood function for bearings-only TMA multiplied by the pdfs (48) for the SNR measurements, after substituting the target parameterization [x n y n ] T = H(t n )¸: (49) The ML estimate is straightforward to compute in this case because propagation loss is due to cylindrical spreading, which has a simple analytical form that makes the gradient of the CAPS model easy to compute. In more realistic situations where the CAPS model does not have a known analytical form (e.g., if it is defined by a set of tabulated values), computing the CAPS model gradient may be difficult. In contrast, the EMAP estimate for this problem does not require computing the gradient of the CAPS model. Fig. 7 shows the estimation errors˜´ˆ¡¸for a 1000 trial Monte Carlo simulation of the EMAP and ML algorithms for this problem. For each trial, new measurements Z = fµ n , » n g N n=1 were generated, and the target parameters were initialized for a target moving from 20 km away at a bearing equal to µ 1 to a location 20 km away at a bearing equal to µ N , measured with respect to the initial and final locations of the sensor, respectively. An increase ² in the logarithm of the joint empirical target pdf of less than 1 £ 10 ¡8 was used as the stopping criterion. The average number of EMAP iterations for each trial was approximately 80, and computing time per iteration was approximately 0.1 s in Matlab on a 600 MHz Pentium III. The 90 percent containment ellipses based on the sample covariances of the EMAP and ML estimation errors are also shown (the ML point estimates are not plotted). Clearly, the ML estimates are well approximated by the EMAP estimates. Fig. 8(a) shows density plots of the missing-range pdfs (31) at several iterations for a typical trial of the EMAP algorithm; thus, each row of the 5 plots comprising Fig. 8(a) sums to one. As the iterations increase, the peaks of these distributions move closer to the target ranges that generated the measurements. Before the first iteration, the peaks are far from the true target ranges because of the (deliberately) poor initialization. After the first iteration, the estimate of¸improves markedly, and the peaks take large steps toward the solution. The last plot in Fig. 8(a) shows that the algorithm converges by iteration 25. Fig. 8(b) shows plots of the missing-range pdfs (31) at the same iterations for a simple down-range variance inflation/deflation strategy whereby the initial value is increased by a factor of 3 for the first 5 iterations, and then decreased to the same value of · used in Fig. 8(a) after the 5th iteration. This simple "annealing" schedule for the down-range variance enables the algorithm to essentially converge after the 5th iteration. This example gives compelling evidence that the EMAP algorithm can easily be made to converge very rapidly. Fig. 9 shows the logarithm of the joint empirical target pdf at each iteration for the example in Fig. 8(a). This plot suggests that the convergence criterion ² = 1 £ 10 ¡8 is overly strict and should be relaxed. This convergence behavior, marked by fast initial convergence, followed by slow steady-state convergence, is typical of EM algorithms in our experience.

VII. CONCLUDING REMARKS
A new formulation of the augmented bearings-only TMA problem is presented. The formulation is novel in that it is based on an empirical MAP approach where measurements are assumed independent a priori because they are predetection measurements; that is, measurements are assumed unconditionally independent until proven conditionally independent by a detection decision. Thus, the EMAP approach is a joint track detection-estimation method.
An algorithm for solving the augmented bearings-only TMA problem using the EMAP formulation is derived by the method of EM. The EMAP target parameter estimation algorithm is an iteratively reweighted, linear least-squares algorithm, provided source level is independent of aspect angle. Thus, the EMAP algorithm is a linear algorithm in most cases of practical interest. It is a remarkable corollary of the EMAP formulation that the classical (i.e., nonaugmented) bearings-only TMA problem can be solved by an iteratively reweighted, linear least-squares algorithm.
The EMAP approach leads naturally to an integral representation of the traditional measurement density in which possible lack of fidelity in the CAPS model is compensated by adjusting the down-range variance · 2 of the geometric kernel. In effect, the proposed compensation averages the CAPS model prediction over a sliding Gaussian window whose size, both down-range and cross-range, increases linearly with range. It is shown that, in the limit, as the down-range length · of the averaging window goes to zero, the new EMAP approach is equivalent to the standard ML approach.

APPENDIX A. ML AS ASYMPTOTIC LIMIT OF EMAP
Laplace's method for asymptotic expansions of integrals is used in this section to reveal the close relationship between the empirical target pdf (18) and the likelihood function (19) of bearings-only TMA. The generic Laplace-type integral is where Á(t) is such that its absolute minimum on the interval [a, b] occurs at the point t = t 0 , where a < t 0 < b, Á 0 (t 0 ) = 0, and Á 00 (t 0 ) > 0: It is assumed that f(t) has at least two and Á(t) has at least four holds as À ! 1. It is derived rigorously using Watson's Lemma in [15]. The second term in the expansion is given in [2]. The particular integrals of interest here take the form (52) = 1, 2, 3, (see (37)), where a(r) = 1 2 µ 1 ¡ x cos µ + y sin µ r ¶ 2 (53) Fig. 9. Logarithm of joint empirical target pdf (eq. (40)) at each iteration corresponding to example in Fig. 8(a). and b(r) = 1 2 µ x sin µ ¡ y cos µ r ¶ 2 : The special case J 1 (·) is equivalent to (4), as is seen by letting r = r n r 1 = r (n) min p(r) = p » n jr n µ n s − n (» n j r n , µ n , s − n ; ³ X n ) and by algebraically manipulating the exponent in (4). The term c´c(Z n , s − n , ³ X n ) is the proportionality constant in (4). The asymptotic form of (52) is sought as · ! 0; hence, · ¡2 plays the role of À, and a(r) plays the role of Á(t) in (50). The necessary condition for the minimum of a(r) is that its derivative be zero. The unique root of a 0 (r) = 0 is r 0 = x cos µ + y sin µ. It is assumed that r 0 is interior to the range of integration in (52). The second-order condition is also satisfied; that is, The function plays the role of f(t) in (50). Using (51) and the fact that a(r 0 ) = 0 gives the asymptotic result J`(·) = g(r 0 ) s 2¼ · ¡2 a 00 (r 0 ) , Substituting (56) and (57) into (58) and simplifying gives J`(·) = % p 2¼¾ p(x cos µ + y sin µ) jx cos µ + y sin µj`¡ 1 , · ! 0: (59) The result (59) for`= 1 is closely related to the pseudolinear approximations used for bearings-only TMA. The earliest references to approximations of this type are [16,17]; later references include [18,19,20]. Pseudolinear methods are discussed extensively in [21].
Because (x cos µ + y sin µ, x sin µ ¡ y cos µ) are the coordinates of the point (x, y) after rotating the coordinate system by µ, x sin µ ¡ y cos µ x cos µ + y sin µ = tan h tan ¡1 ³ y where the approximation in (60) is valid for small bearing measurement errors. Substituting (60) into (59) gives the approximation p(x cos µ + y sin µ) jx cos µ + y sin µj`¡ 1 For`= 1, (61) is The result (62) is important because it shows that the traditional ML formulation of the bearings-only TMA problem is recovered in the limit as · ! 0, after marginalizing over the SNR measurements. Greater accuracy may be sought by using the next term in the asymptotic expansion; however, evaluating the next term requires computing derivatives of the CAPS density. The necessary derivatives may be computed for particular examples, but such derivative computations may be of little practical use in general.
Care must always be exercised in the use of asymptotic expansions to replace numerical integrals, and this application is no exception. Applying the result (59) to the integrals (37) gives the asymptotic ratios for`= 2,3 n ¡ ³ x n ) cosµ n + (y − n ¡ ³ y n ) sin µ n j`¡ 1 , However, if these first-order results are substituted into the EMAP algorithm in place of the numerical integrals, the dependence of the EMAP iterations on the CAPS model is lost entirely. Thus, higher order terms are needed if the asymptotic expansion is to be used in the EMAP algorithm.

APPENDIX B. EMAP AUXILIARY FUNCTION DERIVATION AND MAXIMIZATION
The EMAP auxiliary function for augmented bearings-only TMA of a constant velocity target is derived in this Appendix. General discussions of the EM method are widely available (see [13,22,23]). The joint empirical target pdf (24) is a particular example from a general class of infinite Gaussian mixture densities. For applications of the EM method to finite Gaussian mixtures, see [24][25][26]. The EMAP algorithm may be derived by discretizing the integrals in (24) (see [2]); however, discretization needlessly obscures the discussion. Instead, integrals are retained in the following derivation.
The incomplete-data pdf for this problem is the joint empirical target pdf (24), rewritten here as L inc (Z; ³ X,¸) = N Y n=1 Z r (n) max r (n) min p » n jr n µ n s − n (» n j r n , µ n , S(°n(¸)); ³ X n ) £ N(H(t n )¸j1 n (r n , µ n ), r 2 n § n (µ n ))r n dr n : (64) The missing data in this problem are the range measurements R = fr n g N n=1 associated with the observed bearing and SNR measurements Z = fµ n , » n g N n=1 . The complete data is the union of these two sets; thus, the complete-data pdf is L com (Z, R; ³ X,¸) = N Y n=1 p » n jr n µ n s − n (» n j r n , µ n , S(°n(¸)); ³ X n ) £ N(H(t n )¸j1 n (r n , µ n ), r 2 n § n (µ n )): (66) where the Jacobian J(R) = r 1 ¢¢ ¢r N and dR = dr 1 ¢ ¢¢ dr N , it is easily verified that the marginal pdf of (65) with respect to the missing data R is equivalent to the observed-data pdf (64). Using Bayes Theorem, the missing-data pdf is given by where the missing-data densities p miss(n) (r n j Z n ; ³ X n ,¸) are given by (31). The following two identities are obtained from (68) and (31), Z p miss (R j Z; ³ X,¸)J(R)dR = 1 (69) Z p miss (R j Z; ³ X,¸)J(R n r n )d(R n r n ) = p miss(n) (r n j Z n ; ³ X n ,¸) where the integral over R n r n is similar to (66), but with the integral over r n omitted. The auxiliary function of the EM method is defined by the expectation of the logarithm of the complete-data pdf with respect to the missing data, conditioned on the incomplete data and parameterized by a given estimate for the parameters to be estimated, denoted by¸0. For this problem, the auxiliary function is defined by ª (¸)´Z [log L com (Z, R; ³ X,¸)]p miss (R j Z; ³ X,¸0)J(R)dR: (71) Substituting (65) into (71), interchanging the order of the summation and the integration, and using the identity (70) gives ª (¸) = N X n=1 Z r (n) max r (n) min logfp »njrnµns − n (» n j r n , µ n , S(°n(¸)); ³ X n ) £ N(H(t n )¸j1 n (r n , µ n ), r 2 n § n (µ n ))r n g £ p miss(n) (r n j Z n ; ³ X n ,¸0)r n dr n : Taking the logarithms in (72) and dropping terms in the resulting expression not dependent on¸gives the auxiliary function as the sum of three integrals; see (27)-(30). The M-step maximizes the auxiliary function ª (¸) with respect to¸to obtain the updated estimatȩ + . Each EM step guarantees an increase in the incomplete-data pdf, so successive applications of EM converge to a local maximum if the incomplete-data pdf is bounded above. This monotonicity property is an excellent check against implementation error.
In the important special case when target source level is independent of aspect angle, the last two terms in (27) can be neglected, as they are no longer functions of¸. The M-step for this case becomeş + = argmiņ N X n=1 Z r (n) max r (n) min p miss(n) (r n j Z n ; ³ X n ,¸0)Q(r n , µ n ;¸) dr n r n (73) where definition (26) has been used for the quadratic form Q(r n , µ n ;¸). Differentiating ª MSE (¸) with respect to¸and setting the resulting expression to zero gives the necessary conditions for the minimization in (73), n (µ n )fH(t n )¸+ ¡1 n (r n , µ n )g £ p miss(n) (r n j Z n ; ³ X n ,¸0) dr n r n = 0: Equation (74) satisfies the sufficient condition foŗ + to be a minimum; that is, the matrix of second derivatives of ª MSE (¸) with respect to¸is positive definite.
Equation (74) reduces to the linear system given by (33)-(35). Thus, this linear system holds for all CAPS models when target source level is independent of aspect angle; in this case the CAPS model only contributes to the missing-data densities p miss(n) (r n j Z n ; ³ X n ,¸0).