Misspecification of Multiple Scattering in Scalar Wave Fields and its Impact in Ultrasound Tomography

— In this work, we investigate the localization of targets in the presence of multiple scattering. We focus on the often omitted scenario in which measurement data is affected by multiple scattering, and a simpler model is employed in the estimation. We study the impact of such model mismatch by means of the Misspecified Cram´er-Rao Bound (MCRB). In numerical simulations inspired by tomographic inspection in ultrasound nondestructive testing, the MCRB is shown to correctly describe the estimation variance of localization parameters under misspecification of the wave propagation model. We provide extensive discussion on the utility of the MCRB in the practical task of verifying whether a chosen misspecified model is suitable for localization based on the properties of the maximum likelihood estimator and the nuanced distinction between bias and parameter space differences. Finally, we highlight that careful interpretation is needed whenever employing the classical CRB in the presence of mismatch through numerical examples based on the Born approximation and other simplified propagation models stemming from it.


I. INTRODUCTION
In order to accurately post-process ultrasound measurement data, it is necessary to consider the intricate effects that arise during propagation.Such effects depend on physical parameters of the corresponding measurement scenario.Characteristic examples that have drawn the attention of researchers include shadowing, wherein a large defect obscures the presence of another [1], [2]; reverberation, in which multiple scattering produces repeated echoes that clutter the data and make it difficult to interpret [1], [3]; and attenuation and dispersion, whose impact varies in severity based on the medium characteristics such as grain size in metals [4]- [6].
Among the signal processing community, the phenomenon of multiple scattering has been of particular interest.On one hand, multiple scattering is seen as one of the root causes of image clutter and its negative effects must be alleviated in post-processing [3], [7], [8].On the other hand, some authors posit that the presence of multiple scattering can be beneficial, enabling more accurate estimates of the parameters of interest [9], [10], and going as far as to propose measurement systems in which additional scatterers are introduced artificially in order to promote multiple scattering [11].In both cases, practical data models and estimators are chosen with simplicity and computational efficiency in mind.In cases where multiple scattering is seen negatively, these practical models and estimators are applied to measurement data that follows a different model than the assumed one.On the contrary, when theoretical lower bounds are derived, it is often under the crucial assumption that the chosen data model perfectly matches the true model from which the data was drawn.
More concisely: the two views are at odds due to model mismatch that is often ignored when deriving theoretical performance results.Practical data models are often implicitly or explicitly based on the first order Born approximation, under which multiple scattering is ignored [12].This simplifying assumption results in models and estimation algorithms that are simple and tractable, but it is only valid in limited regimes [13], [14].In practice, these models and estimators are applied to measurement data in order to formulate and solve an inverse problem that yields an image and/or estimates of physical parameters.As the data does not exactly match the model, the interpretability of the images and parameters estimates is negatively affected [15].The applicability of theoretical performance results in these practical, mismatched scenarios depends on whether the mismatch is included in the analysis.
The works [9]- [11], [16], [17] present results on the impact of multiple scattering on the estimation of parameters from scattered wave fields.The authors rely on the Cramér-Rao Bound (CRB) as their tool of choice, as it provides a lower bound on the variance that any unbiased estimator can achieve, and its computation is based only on the data model [18].Theoretical performance bounds of this kind are valuable in experiment design and benchmarking of common systems involving parameter estimation [19]- [23].However, the usage of the classical CRB comes with the assumption that the data model used in the estimation task perfectly matches the true model underlying the data.In practice, this assumption translates into the two following cases.If multiple scattering is considered in the CRB, it is present both in the observed data and the fitted model.If multiple scattering is not considered, it is absent from both.These cases are mutually exclusive, meaning a direct comparison between the two may not provide exploitable knowledge.In addition, neither case considers the more practical scenario in which a model without multiple scattering is chosen for post-processing when multiple scattering is present in the data, in which case the model is misspecified.
One way to alleviate the gap between theoretical performance bounds and the difficulties found when using practical estimators is to entirely replace the models and estimators with ones that are fully physically motivated, removing the mismatch.The field of seismology pioneered this approach through their Full Waveform Inversion (FWI) techniques in which the estimators directly involve the wave equation [24], [25].More recently, medical [26]- [29] and Nondestructive Testing (NDT) [30], [31] ultrasound tomography techniques based on the same principle have been proposed.
Although the FWI framework has the potential to greatly improve imaging quality, it comes at the cost of great computational complexity.FWI cannot be performed analytically, instead relying on iterative algorithms.Accurate estimates are obtained based on approximations to the deviation between the physical model and the observed data.This is necessary due to the complex, nonlinear nature of the model and, consequently, the highly non-convex nature of the target function being optimized.Furthermore, the forward model is computationally demanding.As a result, proper hardware infrastructure is required and reconstruction is time-consuming.What about cases in which such computational effort is inadmissible, e.g. when reconstructions must be done in near real-time or data volumes are overwhelmingly large, or scenarios in which the practice hints at simpler models yielding sufficiently good estimation results?In this work, we employ the Misspecified Cramér-Rao Bound (MCRB) as a tool to analyze the misspecification of multiple scattering.
The MCRB is an extension to the CRB that explicitly accounts for model misspecification.The framework for the MCRB was established more than 50 years ago.However, it is only recently that it has drawn more interest thanks to a wave of landmark publications regarding the MCRB and its usage [32]- [38].Therein, the authors present the same key motivation behind the usage of the MCRB: there are scenarios where, even though the true model may be available, using it directly is prohibitively expensive.This has sparked a variety of works investigating the impact of model misspecification, e.g. in time delay [39] and frequency [40] estimation.Sharing this interest, we turn our focus to the scenario where a misspecified model is chosen on the grounds of being tractable, with the goal of studying the impact of this choice on the achievable performance.
The goal of this work is to study and discuss the impact of the misspecification of multiple scattering on the achievable point scatterer localization performance through the lens of the MCRB, and to do so in a manner that is accessible to the ultrasound NDT and signal processing communities.In the interest of self-containedness, we review key elements of the MCRB and the modeling of ultrasound data.However, the analysis presented in this work is transferrable to other localization tasks based on scalar wave fields, e.g.medical ultrasound imaging and radar localization.
We employ a generalized Slepian formula [32], [36], [38] for the MCRB of misspecified-unbiased estimators.We focus on the estimation scenarios where the assumed model has a misspecified mean µ : R d → C N which differs from the true mean s ∈ C N , and both are afflicted by circularly symmetric complex Gaussian noise with a known covariance matrix Σ ∈ C N ×N .In particular, the true data model follows the scalar Helmholtz equation, while the misspecified model is the classical time delay model based on the first order Born approximation, consisting of scaled and shifted copies of a predetermined pulse shape.We then perform numerical experiments for an ultrasound NDT computed tomography scenario consisting of two point-like scatterers, noting that the results obtained are valid for any number of scatterers so long as they are only pairwise close to each other.Finally, we compare the MCRB to the classical CRB and discuss the implications of the true and misspecified models having different parameter spaces, limiting the usefulness and interpretability of the classical CRB.

II. MISSPECIFIED CRAMÉR-RAO BOUND
We begin by briefly discussing the MCRB and establishing the key components needed for its computation.For in-depth derivations of the MCRB and discussions on the class of estimators to which it applies, we refer the reader to [32]- [38], [41] and the references therein.The definitions and formulations presented next are based on the aforementioned references.
A. GENERAL DEFINITIONS Denote by x ∈ C N a single observation obeying an unknown probability density function (PDF) p X (x), shortened as p X and referred to as the true PDF.Instead of this true PDF, the observation x is assumed to follow a PDF f X (x|θ), shortened as f X |θ , which depends on the parameter vector θ ∈ Θ ⊂ R d .Throughout this work, we refer to the PDFs as statistical models.If there is no θ for which p X (x) = f X (x|θ), the assumed model is said to be misspecified.These conditions set the context for the following definitions and results.
Definition II.1 (Pseudo-True Parameter).Consider the Kullback-Leibler Divergence (KLD) from f X |θ to p X defined as If there exists a unique parameter θ 0 ∈ R d such that then θ 0 is referred to as the Pseudo-True Parameter (PTP).
The PTP in Definition II.1 can be interpreted as the parameter of the misspecified model that best explains data stemming from the true model.We are interested in estimators that relate to the PTP in the following way.

Definition II.2 (Misspecified Unbiasedness
).An estimator θ of θ 0 based on the misspecified model f X |θ is said to be unbiased in the misspecified case, or Misspecified Unbiased (MS-unbiased), if and only if where E p represents expectation with respect to the true PDF p X .
The goal is to obtain a lower bound for the error covariance of MS-unbiased estimators following Definition II.2.To this end, the following auxiliary quantities are necessary.
Based on the MLL, the auxiliary matrices defined entrywise as and are of special interest.In the case of correctly specified models, the matrices defined in ( 5) and ( 6) are equivalent to each other (up to their sign) and correspond to the Fisher Information Matrix (FIM).
The MCRB can now be defined based on the preceding definitions by noting that, in the misspecified case, matrices (5) and ( 6) are no longer equivalent.Instead, their non-equality is central to the MCRB.
Definition II.4 (Misspecified Cramér-Rao Bound).Let f X (x|θ) be misspecified with respect to p X (x).Let θ be an MSunbiased estimator of θ 0 following Definition II.2.Consider additionally the error covariance matrix of θ given by If the matrix A θ0 defined entrywise as in ( 5) is invertible, then where a matrix inequality of the form U⪰V signifies that U − V is positive semidefinite and MCRB(θ 0 ) is the Misspecified Cramér-Rao Bound.
To compute the MCRB, a choice of true and misspecified models must be made and used in tandem with these definitions.Note that (8) depends explicitly on the PTP.In general, computing the PTP analytically is challenging.The following theorem allows this to be circumvented.

Theorem 1 (Mismatched Maximum Likelihood Estimator).
The quantity is referred to as the Mismatched Maximum Likelihood Estimator (MMLE).If regularity conditions are met, the MMLE is asymptotically (as the number of observations goes to infinity) MS-unbiased and its error covariance coincides with the MCRB.Proof: The works [33], [34] discuss this result in detail and refer to the original proof presented by Huber and White.■ Based on Theorem 1, the MMLE is a practical alternative to the analytic computation of the PTP and is the method of choice in this work.We furthermore narrow our focus to the study of Gaussian models with misspecified mean, for which specialized formulations of the KLD, MLL, and MCRB can be obtained.These formulations are presented next.

B. GAUSSIAN MODELS WITH MISSPECIFIED MEAN
Let y ∈ C N be the observed data which is corrupted by zero mean, circularly symmetric, additive complex Gaussian noise n ∈ C N ∼ CN (0, Σ), where Σ ∈ C N ×N is the noise covariance matrix.The observed data can then be written as y = s + n.The vector s ∈ C N is the true mean of the data, and so y ∼ CN (s, Σ) defines the true PDF p Y (y).Misspecification is introduced through the assumption (be it because s is unknown or due to practical considerations) that the data y has a mean µ(θ) which is different from s for all values of the parameter θ.This means that the data is assumed to obey y ∼ CN (µ(θ), Σ), which in turn defines the misspecified PDF f Y (y|θ).
In this scenario, the KLD from f Y |θ to p Y takes the wellknown form Similarly, the MLL is given by v(θ) = − ln(π N det{Σ}) Applying the definitions from Section II-A to this Gaussian model scenario, a generalization of the Slepian formula [18] is obtained.We employ the formulation shown in [32], dubbed Deterministic Generalized Slepian Formula in [38], noting that the resulting bound applies to any MS-unbiased estimator [33]- [35], [41].The expressions are reworked for the specific case where the model parameters are real, i.e. the mean is a function µ : R d → C N .
Theorem 2 (Generalized Slepian Formulae).Assume the regularity conditions hold so that the MCRB is given by MCRB(θ 0 ) . Consider a Gaussian model whose true mean s ∈ C N is misspecified as µ(θ) ∈ C N with θ ∈ R d .The matrices A θ0 and B θ0 are defined entrywise by the generalized Slepian formulae and where Σ ∈ C N ×N is the noise covariance matrix shared by the true and misspecified models and ∆(θ) = s − µ(θ).
Proof: Different derivations are presented in [32], [36], and [38].■ Note that, if µ( θ) = s for some value θ of θ, the model is no longer misspecified and all terms involving the difference ∆(θ) = s − µ(θ) become 0. In this case, B θ0 = −A θ0 , both Slepian formulae correspond to the correctly specified FIM, and the MCRB is equal to the classical CRB.
As the goal is to apply these results to true and misspecified models related to scalar wave phenomena, we discuss the basics of scalar wave fields and wave propagation models next.In the remainder of this work, the upcoming propagation models will be referred to simply as models, as they describe how the means µ and s of the misspeficied and true statistical models, respectively, are constructed.

III. SCALAR WAVE FIELD MODELS
Wave phenomena can be modeled through the Helmholtz equation [12], [14], [26], which in the presence of sources takes the inhomogeneous form ∇ 2 ψ t (ω, x) + k 2 (ω, x)ψ t (ω, x) = −s(ω, x).(14) In ( 14), x ∈ R 2 denotes a position vector, ω is the angular frequency, ψ t : R × R 2 → C is the total field in a region of interest (ROI), k : R × R 2 → C is the wavenumber, and s : R×R 2 → C is a harmonic excitation source.The wavenumber can be decomposed into where c : R 2 → R assigns a propagation speed to each point in space and β : R×R 2 → R is a frequency and space dependent attenuation factor.
It is important to note that the properties of β depend on the application: when dealing with ultrasound and human tissue, β is approximately linearly dependent on the angular frequency ω [4]; however, a nonlinear dependence on ω may be present when dealing with ultrasound in metals [5], [6].Assuming a linear dependence on ω for simplicity, the wavenumber k can be rewritten as with an attenuation factor β that depends only on the position x after factoring ω.Problems involving (14) are categorized into forward problems where k(x) and s(ω, x) are known and ψ t (ω, x) is sought after, and inverse problems in which k(x) is estimated given s(ω, x) and samples of ψ t (ω, x).A common solution approach for inverse problems is to introduce a contrast term so that only deviations with respect to a reference value of k are considered [12], [14], [26], [27], [42], [43].One possible formulation of the contrast is where γ : R × R 2 → C is referred to as the wavenumber contrast, and k 0 ∈ C is a complex scalar denoting the background wavenumber from which k deviates.By recalling the assumption that k depends linearly on ω, the quantity can be defined, allowing the contrast γ to be reformulated as Due to this relationship, the quantity γ : R 2 → C as given by ( 18) is referred to as the frequency flat contrast.The frequency flat contrast allows ( 14) to be rewritten in the form (20) where ∇ 2 and k 2 0 have been grouped together to highlight their behavior as a differential operator acting on ψ t .Crucially, this A transducer ring is placed inside a medium that extends infinitely in all directions.Within the medium, a darker gray color represents a wavenumber that differs from that of the background.The total field within the medium can be expressed in terms of two distinct components.The incident field is shown in Figure 1b, where it is clear that the inhomogeneity is no longer present and the excitation produced by a transducer propagates in the homogeneous medium.In Figure 1c, the inhomogeneities are shown behaving as sources producing a field that also propagates through the homogeneous background medium characterized by k0.
operator now depends on k 0 rather than the spatially varying k.Paired with the Sommerfeld radiation condition, this allows (20) to be solved by using the Green's function associated to the differential equation.The Green's function where H (2) 0 is the zeroth order Hankel function of the second kind, and x, x ′ ∈ R 2 denote locations in the ROI.
Using the Green's function in (21), the total field can be expressed as with subindices i and s referring to incident and scattered fields, respectively, and where dependence on ω has been omitted.In turn, these fields are given by and where Ω refers to the support of an excitation source or transmitter and Γ is the ROI.As (23) does not depend on γ, the incident field ψ i is considered to be known or computed ahead of time.Equations ( 23) and ( 24) can be interpreted as decomposing the total field into the contributions from two kinds of sources: the original sources that excite a medium where γ is everywhere zero, and regions in the medium where γ ̸ = 0 which behave as additional sources.This is illustrated in Figure 1.Note that, since the contrast γ behaves as a source in (24), the scattered field scales with the magnitude of the contrast.Equations ( 22)-( 24) can be discretized through different choices of sampling grids and integral approximations.These choices result in distinct forward models, as discussed next.

A. HELMHOLTZ MODEL
Consider an ROI sampled with resolution ∆ x along the horizontal axis and ∆ z along the vertical axis so as to yield N z • N x points on a regular 2D grid.After discretization, we represent the total and incident fields in the ROI at frequency ω and the frequency flat contrast as vectors of the form ψ t , ψ i , γ ∈ C Nz•Nx .After choosing a numerical integration scheme for (24), a matrix G Γ Γ ∈ C Nz•Nx ×Nz•Nx is defined so that the total field in the ROI defined in ( 22) can be approximated as where the factor ω 2 has been absorbed into G Γ Γ , • denotes the Hadamard or elementwise product, and diag{•} yields a diagonal matrix whose entries correspond to the given argument.Furthermore, in the absence of self-interference and considering point-like receiver and transmitter elements, the field observed at a receiver is equivalent to the scattered field given by ( 24) when evaluated at the coordinates of said receiving element.For an array with M R receiving elements, a matrix G A Γ ∈ C MR×Nz•Nx can be defined so that ( 24) is approximated by by evoking the commutativity of the Hadamard product to swap ψ t and γ, and where ψ s ∈ C MR is the scattered field measured at the receiver array at frequency ω.Data comprising N ω frequency bins of interest and M T transmitting elements can be obtained by using ( 25) and ( 26) N ω • M T times.We refer to ( 25) and ( 26) as the Helmholtz model throughout this work.Note that, in the special case where G Γ Γ and G A Γ are constructed by simply taking samples of the Green's function, (25) and ( 26) correspond to the Foldy-Lax model [9]- [11], [44], [45].This is equivalent to considering the contrast γ as being caused by point scatterers represented by a sum of delta functions in space.Substituting such a contrast term into (22) and ( 24) directly yields ( 25) and ( 26).

B. BORN MODEL
In both forward and inverse problems involving the Helmholtz equation ( 14), the largest computational effort is spent in solving (25) for ψ t .A common solution is to make the simplifying assumption that, at the scatterers, the total field on the right-hand side of ( 25) is equivalent to the incident field.The justification behind this is that, if the contrast γ is small, then the incident field dominates the scattered field and ψ t = ψ i + ψ s ≈ ψ i .This is known as the first order Born approximation, here referred to as the Born model, and it can be interpreted as ignoring the interactions among the scatterers and therefore omitting multiple scattering.A full data set can then be obtained by computing for all combinations of frequency and transmitter indices.

C. DELAY MODEL
A further simplification can be made by considering point-like scatterers and modifying the Green's function in two ways.First, the Green's function can be approximated as for large distance values ∥d∥ 2 = ∥x − x ′ ∥ 2 .Next, the k 0 ∥d∥ 2 in the denominator is dropped.This amounts to ignoring the beam spread as waves travel away from their source.If we assume the background wavenumber k 0 to be real, we assume the medium to be non-attenuating.This turns (28) into the transfer function of a time delay.The overall forward model subsequently turns into a sum of copies of a transmitted waveform, each of which has been delayed depending on the geometry concerning the sensors and the scatterers, and each copy is additionally scaled based on the value of γ.
These simplifications can be used to rewrite the forward model as follows.As the transmitters are considered to be point-like, the excitation s(ω, x) is of the form a(ω)δ(x − x t ) for the t th transducer located at x t .The amplitude a : R → C corresponds to the spectral amplitude or Fourier coefficient at angular frequency ω of an excitation pulse shape p(t) with spectrum P(ω).Scaled, delayed copies of this pulse shape are measured at the receivers, one for each scatterer.The scattered field observed at a fixed receiver after a medium with U scatterers is excited by a transmitter is then approximately given by where q u ∈ C is the scattering coefficient of the u th scatterer and τ (x u ) is its associated time delay, noting that this time delay depends both on the scatterer's location x u and the location of the transmitter and receiver.The time delay can be written as for a transmitter located at x t and a receiver at x r .Such a model is employed in migration and delay and sum techniques for imaging and localization, e.g. the Synthetic Aperture Focusing Technique [46] and the Total Focusing Method [47] in ultrasound.A comparison among the three aforementioned models is shown in Figure 2 in the time domain.In this figure, the Helmholtz model refers to fields produced via ( 25) and ( 26), the Born model corresponds to (27), and the Delay model refers to (29).The plot shows the scattered field observed by a point receiver when a medium containing two point scatterers is excited by a point transmitter emitting a Gabor function.The exact details of the simulation are presented in Section V. Note that the Helmholtz model, shown in black, exhibits multiple scattering which manifests as reverberation.The Born approximation, in blue, agrees with the two largest pulses that can be observed in the Helmholtz model; however, reverberation is absent.Finally, the delay model lacks reverberation and, although the amplitudes of the largest echoes are similar to those of the previous two models, there is also a discrepancy due to the simplified Green's function.

IV. MISSPECIFICATION OF MULTIPLE SCATTERING
We now apply the MCRB from Section II to study the achievable estimation performance in localization tasks involving the scalar wave phenomena from Section III.Misspecification occurs when estimating parameters based on a simplified model such as the Born or delay model, ignoring the presence of multiple scattering in the data.
We focus on the scenario in which the model follows a circularly symmetric complex Gaussian distribution whose mean is misspecified.We highlight the interpretability of the generalized Slepian formulae presented in Section II-B, which explicitly account for the mismatch between the true and assumed models.More concretely, when applying the generalized Slepian formulae in the present application, the true and misspecified means correspond to different scalar wave propagation models, and the MCRB is explicitly formulated in terms of their difference.

A. MISSPECIFICATION OF THE PROPAGATION MODEL
We contribute to the discussion on the impact of multiple scattering on estimation performance by analyzing the model mismatch via the MCRB.An MCRB that evaluates the impact of the misspecification of multiple scattering can be obtained by combining the models in Section III and the generalized Slepian formulae in Theorem 2. In order to account for multiple scattering, the true mean is taken to be s = ψ s ∈ C Nω•MR•MT as given by the solution to the Helmholtz equation obtained by solving equations ( 25) and ( 26) for all N ω frequencies and M T transmitters.
This true model depends on γ, which has the locations x u of U point-like scatterers as its support and assigns the scatterers a wavenumber contrast γ u = γ(x u ) with respect to the background medium.For the sake of clarity, let us define a true parameter ϕ containing the aforementioned scatterer locations x u and their corresponding wavenumber contrasts γ(x u ).We can then write the true mean as s(ϕ).The model is additionally corrupted by circularly symmetric complex Gaussian noise n ∼ CN (0, Σ).
The misspecified mean is chosen as the delay-based model from (29) after discretization in the frequency domain, and the same noise n as before is considered.This results in a model µ(θ) = ψ D ∈ C Nω•MR•MT .The parameter θ of the misspecified model describes the locations x u of the scatterers just as ϕ does, but instead of wavenumber contrasts, a complex scattering coefficient q u = a u exp(ȷθ u ) characterized by a realvalued amplitude a u and phase θ u is assigned to each scatterer as in (29).
Crucially, the parameter spaces of the true and misspecified models are overall distinct and intersect only in the subspace of the location parameters.In the context of model misspecification, even if the parameter ϕ of the true model s is known, the assumed model µ cannot in general be evaluated at ϕ, e.g. when ϕ and θ are of different dimensions.Even in the case when µ(ϕ) is properly defined, it has no connection to s(ϕ) as the physical meaning of the parameter spaces is different.The scattering coefficients q u in θ represent the combined effect of the corresponding wavenumber contrast γ(x u ), the Green's function, and beam spread, and attribute them to a single scalar acting on the incident field.This is a large simplification of the true physical phenomenon.As an added consequence, the locations x u of the scatterers according to the assumed model µ generally differ from those of the true model.
In addition to the previous observation, recall that the MCRB provides a lower bound for the error covariance matrix C p ( θ, θ 0 ).That is to say, the MCRB describes the variance of MS-unbiased estimators around the PTP θ 0 .The MCRB does not, however, compare the estimated parameter θ against the true parameter ϕ.Such a comparison is difficult, since the quantities ϕ−θ 0 and ϕ− θ, as well as the Mean Squared Error (MSE), can only be well-defined on intersecting subspaces of the parameter spaces.Furthermore, there is no general framework with which the difference between the true and pseudo-true parameters can be studied [33], [34].In the special cases where the parameter difference vector r = ϕ − θ 0 can be defined, as is the case in our application, a connection can be drawn between the MCRB and the MSE as discussed next.
Throughout the remainder of this work, the scenario with U = 2 point scatterers is considered.
T , with (•) s as the projection operator onto the subspace of location parameters, the difference vector r = ϕ s − θ s 0 between the true and pseudo-true parameters can be properly defined.The MSE of an MS-unbiased estimator θ with PTP θ 0 and true parameter ϕ on the subspace denoted by (•) s can then be written as [34] MSE( θs Due to ( 8) and ( 31), it follows that We highlight that, in the context of misspecification, ( 31) is not to be understood as a bias-covariance decomposition.Instead, since the error covariance is defined with respect to the PTP as in (7), a proper bias-covariance decomposition pertains to the Misspecified Mean squared Error (MMSE) and is given by where θ = E p { θ} and with b = θ 0 − θ as the bias.When θ = θ 0 , the estimator is MS-unbiased.As a consequence, the bias b is zero and the MMSE and error covariance C p ( θ, θ 0 ) of the estimator θ coincide.

B. EMPIRICAL VALIDATION OF MODEL USEFULNESS AND BOUNDS
An important observation is that equations ( 8) and ( 33) can be employed both to validate the computed MCRB and to verify the MS-unbiasedness of the implemented MMLE.Additionally, (32) is useful in corroborating the usefulness of the delay model in localization tasks through the lens of the MCRB.Note that the delay-based model and other similar models considering only few high amplitude scattering events are wellestablished and have been studied extensively in applications such as direction of arrival estimation, beamforming, and other fields involving the scattering of wave-like phenomena [44], [45], [48]- [50].As such, this verification through the MCRB is to be understood as another tool in the proverbial toolbox, with the important caveat that the true parameters of the true model must be known.This is only possible in simulations and during calibration.
Under the assumption that the requisite regularity conditions hold, the KLD in (10) can minimized numerically in order to obtain the PTP θ 0 , which can then be substituted into the Slepian formulae to compute the MCRB.The computed MCRB can then be compared against the Empirical Misspecified Mean Squared Error (EMMSE) obtained from Monte Carlo simulations in the presence of noise.The EMMSE estimates the MMSE from (33) and is given by (34) based on N estimates θn of the PTP θ 0 .The empirical mean θ can be defined in a similar fashion.This EMMSE can then be decomposed as in (33) to corroborate that the bias with respect to the PTP θ 0 is negligible.This comparison simultaneously shows whether the chosen estimator is a MMLE and whether the EMMSE follows the MCRB.
An Empirical Mean Squared Error (EMSE) can be defined analogously to (34) and used in (32).This provides an additional means of validating the computed MCRB and chosen estimator: if the estimator is a MMLE and the MCRB is computed correctly, equality is achieved in (32) in the limit as the number of realizations goes to infinity.As an added benefit, the computation of the difference vector r offers insight into the usefulness of the chosen misspecified model.If the entries of r are small, the delay model correctly conveys scatterer location information in spite of ignoring beam spread and multiple scattering.After this twofold validation, the MCRB can be used in further numerical simulations in which noise is accounted for through the noise covariance Σ instead of Monte Carlo trials, reducing the computational load.

V. SIMULATIONS
In the numerical simulations presented next, inspiration is taken from ultrasound NDT.In particular, the simulations concern themselves with the task of locating defects at which the Speed of Sound (SoS) deviates from a known background value.The background medium is set to have an SoS of c 0 = 6400 m s −1 and no attenuation, i.e. β 0 = 0.A total of U = 2 point scatterers are placed in the medium, meaning the frequency flat contrast γ is nonzero at exactly two locations x u as discussed previously.Both defects are set to have the same wavenumber k.
For the excitation s(ω, x) = a(ω)δ(x − x t ) in ( 14) and P(ω) in ( 29), a modulated Gaussian of the form is chosen.The parameters for the pulse shape described in (35) are chosen so that they correspond to a reference real world measurement.The sampling frequency is set to f s = 40 MHz, the bandwidth is controlled by the factor α = (4.67MHz) The geometric parameters are chosen in terms of the reference wavelength λ 0 = c 0 /f c as follows.The scatterers are located within a ROI of size N z × N x = 81 × 81 extending over a square region of size (10λ 0 ) 2 with resolution ∆ z = ∆ x = λ 0 /8.A total of 32 transducers are placed around the ROI along a circle of radius 10λ 0 at regular intervals.An illustration of the geometry is shown in Figure 3.For the sake of simplicity, the transducers do not interact with the medium and therefore produce no further reflections.
The simulation scenarios are constructed by varying the parameters as follows.One defect is moved to each of the possible 81 × 81 positions.The second scatterer is kept fixed at the origin, i.e. x 2 = [0, 0] T .For each of the resulting configurations, the MCRB is computed.In this computation, the true parameters are given by ϕ T , where q u = a u exp(ȷθ u ) is the scattering coefficient of each scatterer and x u = [x u , z u ] T are the scatterer locations.This means that MCRB(θ 0 ) is of size 8 × 8.
In this work, our interest lies in defect localization.In light of this, the quantity is taken as a localization performance indicator.The Q matrix is to be substituted by the classical CRB or any of the quantities involved in ( 8), (32), or (33) (e.g.MCRB(θ 0 ) or rr T ).This means that Q is real-valued and is of size 8 × 8 if it is defined over θ or ϕ, or of size 4×4 if the projection (•) s was employed as is the case with rr T .Furthermore, Q ii with i ∈ {1, 2, 3, 4} refers to the first four diagonal entries of the Q matrix, which in all cases corresponds to the location parameters x 1 , z 1 , x 2 , z 2 of the two scatterers.With this definition, the scalar Q represents a form of absolute error.Specifically, depending on the chosen matrix Q, the Q-value represents a sum of standard deviations, a sum of absolute differences, or a sum of absolute biases, and is equivalent to computing elementwise square roots followed by the traces of ( 8), (32), and (33) over the subspace (•) s of location parameters.This enables the creation of an image of 81 × 81 pixels, referred to as a Q-image, where the location of each pixel corresponds to the coordinates x 1 of the moving scatterer as given in ϕ and the pixel intensity is given by the corresponding Q-value which has units [µm].The individual Q-values will be referred to by the quantity they represent, e.g.''Bias'' meaning that Q = bb T is substituted into (36).

A. COMPUTATIONAL DETAILS
When generating fields according to the Foldy-Lax model, equations ( 25) and ( 26) are treated differently.If G Γ Γ in ( 25) is computed simply by taking samples of the Green's function, the condition number is poor and inversion, including direct inversion, is challenging.Instead, the matrix is constructed following the 2D trapezoid rule.This can be interpreted either as a smoothing of the Green's function, or as letting each defect consist of a group of four closely spaced defects.As no inversion is involved in (26), G A Γ is constructed from samples of the Green's function.
Both the PTP θ 0 and estimates θ are required for the computation of the MCRB.Due to the noise statistics, the minimization of the KLD (10), which yields the PTP θ 0 , takes the nonlinear least squares form Asymptotically MS-unbiased estimates can be obtained by minimizing the negative of the MLL (11), which is equivalent to A two step procedure is adopted in solving for θ 0 and θMML .Note, however, that the procedure described next is not suitable for reconstruction and localization in practice, but is only employed given the amount of prior information available when generating simulations.The pulse shape, the grid where the defects are located, and the number of defects are all available in simulations, allowing the construction of a linear model.A dictionary matrix A ∈ C Nω•MR•MT ×2 is built based on the discretized delay model given all the prior knowledge, letting each column correspond to the signal that would be observed if only a single scatterer were present.The only unknown quantities are the scattering amplitudes q u , which can be obtained as q LS = argmin q ∥y − Aq∥ 2  2 = A † y, (39) where (•) † denotes the pseudoinverse.The entries of q LS ∈ C 2 correspond to estimates of the scattering amplitudes q u .When both the true and assumed model correspond to the delay model, this step suffices to obtain all the model parameters.However, the scatterer coordinates in the delay model in general do not coincide with those of the Helmholtz model.
In the second step, the estimates are refined by using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm [51].The target functions are the nonlinear least squares KLD and MML expressions in (37) and (38).The true scatterer coordinates x 1 , x 2 in ϕ, as well as q LS , are used as initial guesses for the parameters θ.We emphasize that this is possible in simulations, but the true parameters ϕ of s are unknown in practice.The derivatives involved in the BFGS algorithm are computed through automatic differentiation using the JAX library [52].

B. VALIDATION
Two scenarios are considered in the validation of the computed MCRBs and the related Q-images.In the first scenario, the defects are given an SoS c = 1500 m s −1 , and the MMLE is computed following the aforementioned procedure.This is done for 200 noise realizations, after which the EMMSE of the estimates of the localization parameters is used to compute a cross section of a Q-image.Next, the PTP is computed by minimizing the KLD.Using Theorem 2, the MCRB is computed based on the PTP and a noise covariance matrix Σ = σ 2 I, with σ 2 = 3.From this result, the corresponding Q-image cross section is computed.Finally, the PTP is taken as the input parameter of the delay model, which is now treated as the true model.This yields one last Q-image cross section corresponding to the classical CRB to be used as a reference.This procedure is then repeated for c = 5000 m s −1 for 500 noise realizations.The cross sections of the Q-images are constructed by allowing x 1 to vary and keeping z 1 = −0.704mm fixed.The cross sections are shown in Figure 4.
The EMMSE-based Q-values are observed to follow those of the MCRB, advocating for its validity.The empirical bias-covariance decomposition of the EMMSE in Figures 4a  and 4b illustrates that the empirical bias is negligible, meaning the estimator is practically MS-unbiased.Figure 4c shows better agreement with the MCRB than Figure 4d, since the magnitude of the frequency flat contrast γ increases as the difference between the background and defect SoS increases.This increases the scattering amplitude and Signal to Noise Ratio (SNR) relative to the simulation scenarios with a lower SoS contrast.The EMMSE agrees with the MCRB instead of the classical CRB for both SoS values, as is to be expected in the presence of model mismatch.Reiterating, Figure 4 is  4a and 4b show the empirical bias-covariance decomposition of the EMMSE in logarithmic scale.The empirical bias is an order of magnitude smaller than the covariance.Since the EMMSE Q-value is obtained through the Pythagorean addition of the empirical bias and covariance, the influence of the bias is negligible and the covariance coincides with the EMMSE.The EMMSE and covariance additionally coincide with the MCRB.In Figures 4c and 4d, the agreement between EMMSE and MCRB is illustrated in linear scale.Even though more noise realizations are used in the computation of Figure 4d, the EMMSE appears noisier than in Figure 4c due to the lower SoS contrast and consequently lower SNR.The classical CRB is shown as a reference, and the EMMSE is observed to follow the MCRB instead of the CRB.
an application of ( 8) and (33) illustrating the validity of the computed MCRB and the MS-unbiasedness of the estimator.Using the same parameters as in Figure 4, the EMSE and empirical parameter error r over the subspace of defect location parameters are computed next.Figure 5 shows a good match between abs(rr T ) and abs(MSE − MCRB), which is possible for MS-unbiased estimators whose error covariance coincides with the MCRB.This corroborates the findings illustrated in Figure 4. Additionally, the Q-values are noticeably small.Consider the worst case scenario in which the parameter error plotted in Figure 5 is attributed to a single location parameter out of the total four, e.g.all of the error is due to x 1 .Even in this case, the parameter error is two to three orders of magnitude smaller than the value of the parameter, and is therefore negligible except for very closely spaced scatterers.This serves as an alternative corroboration of the widespread success of the delay-based model.As a final remark regarding validation, notice that the location parameter difference in Figure 5b is smaller than that in Figure 5a.This coincides with the notion that the degree of misspecification of multiple scattering is lower when the SoS contrast is smaller, otherwise referred to as the weak scattering regime.
The remaining simulations are carried out in the absence of noise, recalling that the MCRB depends on the noise only through the covariance matrix Σ which can be introduced ex post facto.This allows the construction of the full Q-images for all SoS values.

C. SIMULATION RESULTS
Having validated the numerical computation of the MCRB, we move on to simulations without Monte Carlo noise realizations.We consider the following scenarios.The defects are allowed to have wavenumbers defined by c ∈ {1500 m s −1 , 2700 m s −1 , 3800 m s −1 , 5000 m s −1 }, β = 0. Recalling that the background SoS is c 0 = 6400 m s −1 , scenarios ranging from high to low SoS contrasts, and therefore Helmholtz, c = 1500 0.17  wavenumber contrasts, are considered.As the wavenumber contrast decreases, the magnitude of the scattered field and the effect of multiple scattering decreases, the three models discussed in Section III become more similar, and the degree of misspecification decreases.In order to study this behavior, two true models are considered: the Helmholtz model and the delay model.As was done during the validation, the PTP is taken as the input parameter of the delay model.More details on this choice are presented in Section VI.This yields a total of 8 different simulation scenarios: for each of the two true models (the Helmholtz and delay model), each of the four defect speeds of sound is considered and a corresponding Q-image is computed.
To reiterate, the simulation scenarios consist of all possible combinations of four defect SoS values c ∈ {1500 m s −1 , 2700 m s −1 , 3800 m s −1 , 5000 m s −1 } paired with the delay and Helmholtz models, yielding a total of eight scenarios.In all scenarios, the assumed model during estimation is the delay model, meaning that there is misspecification when the true model is the Helmholtz model.The Q-images for the eight simulation scenarios involving different models and scatterer parameters are presented in Figure 6.The axes are shown in wavelengths, i.e.L/λ 0 for a length L, so as to facilitate the interpretation of the distance between the scatterers.
For a fixed SoS c, the Q-images appear to have an overall similar appearance regardless of the true model.However, the MCRB values are overall lower than those of the CRB.This means that, in the simulation scenarios, the presence of multiple scattering in the data results in better estimation performance in terms of the error covariance of θ, across all values of c and ϕ, than would be predicted by the classical CRB.This observation is to be interpreted carefully, and is discussed in detail next.

VI. DISCUSSION
The simulation scenarios in this work hint at an overall performance improvement in localization accuracy in the presence of multiple scattering, even when the wavenumber contrast is low.At first glance, this appears to contradict findings in the literature which state that the presence of multiple scattering may be beneficial or detrimental depending on the configuration of the sensors and scatterers [9]- [11], [16], [17].However, these works consider single transmitter and receiver scenario with harmonic excitation, motivated by the search of analytic expressions.In contrast, our simulations employ a sensor array and a realistic signal with a broad bandwidth, and we study the impact of model mismatch directly.Although we don't present analytic results, the overall procedure for the computation of numerical results and their interpretation can be applied to any scenario.
More importantly, when mismatch is considered in these works, it is studied through the MSE.This automatically introduces a parameter difference r that, in these works, is referred to as a bias.As discussed in Section IV-A, this difference should not be interpreted as a bias in the context of misspecification.Instead, the presence of misspecification turns the MLE into the MMLE.In the context of misspecification and the MCRB, the results reported in [9]- [11], [16], [17] contain a combination of parameter space differences and bias with respect to the PTP.The interpretation of the findings in the present work and a comparison to the aforementioned prior works constitutes the remainder of this discussion.
Before jumping into the discussion, however, the question arises as to why the CRB is consistently higher than the MCRB.Although additional factors may be involved, the connection between total signal energy and multiple scattering can be highlighted.In the present investigation on the misspecification of multiple scattering, model mismatch directly affects the SNR.As was illustrated in Figure 2, the delay model often matches the first scattering event in the Helmholtz model.However, further scattering events are ignored, meaning that the energy content of signals following these models is different.Normalization and scaling are possible, but this difference in signal energy is inherent to the model choice.The presence of multiple scattering can thus increase the SNR when distinct echoes do not overlap or when they interfere constructively with each other.With this observation out of the way, we first address the matter of parameter differences and bias.

A. BIAS OR PARAMETER DIFFERENCE?
Moving onto the interpretation of the results, we reiterate the accuracy of the MCRB when describing the MMSE, and not the MSE.This observation is crucial because, in the presence of misspecification, the MSE involves the computation of differences on intersecting subspaces of the parameter spaces of distinct models.This motivates the usage of the projection (•) s onto the subspace of location parameters when considering the MSE throughout this work.Instead, the MCRB as defined in Definition II.4 incorporates misspecification by considering the error covariance with respect to the PTP, i.e. within a single parameter space.
In the delay model, only the scatterer location parameters overlap with the parameter space of the true model.This means that, for a true model s Helmholtz (where the subscript has be added for clarity) with true parameter ϕ, the quantity µ delay (ϕ) has no useful interpretation.In contrast, the aforementioned works use the Born approximation, which shares its parameter space with the Helmholtz model.It may then appear more reasonable to evaluate the Born model at the true parameter ϕ, but the same property still holds: even though µ Born (ϕ) can be properly evaluated and appears similar to s Helmholtz (ϕ), the two are not directly related since the models are distinct.This can be observed e.g. in [17], where it was noted that the usage of the Born model during estimation introduces a ''bias''.The observed bias was computed based on the MSE and therefore actually refers to the parameter difference r, which in their case is properly defined over the entire parameter space without needing the projection (•) s .The observed parameter difference means that, although the models share the same parameter space, the presence of misspecification in the model has affected the estimation procedure.As mentioned previously, this can be interpreted as misspecification ''turning the MLE into the MMLE'', which we make precise as follows.
Recall that the Mismatched Maximum Likelihood Estimator under misspecification of multiple scattering was formulated in (38) This estimator is asymptotically MS-unbiased with expected value θ 0 .What happens if we now change the model s to the delay model without modifying the noise statistics?Doing this yields the expression in which θ * is the true parameter of the new true model ψ D .This expression clearly has no misspecification and describes the Maximum Likelihood Estimator.When only the true mean of the data is changed, the MMLE turns into the MLE.As such, we can say that ''In the presence of misspecification, since the MLE turns into the MMLE (in the sense described previously), the estimation task automatically deals with the search of the PTP θ 0 , and not the search of the true parameter ϕ of the true model.In ultrasound localization, the presence of model misspecification means that the estimated locations and scattering amplitudes generally do not match the ground truth even if the misspecified model appears to have the same parameter space as the true one.''This directly explains the observed ''bias'', which as mentioned previously is technically the parameter difference r introduced by the misspecification.We can then add that ''The magnitude of the resulting parameter difference r (when it is defined) determines whether the chosen misspecified model is useful for the estimation task.In the context of ultrasound localization, a small parameter difference means defects can be located accurately despite the model mismatch.''We follow up on this point by studying additional properties of the MLE and MMLE.

B. USING THE PTP AS A REFERENCE
Next, the asymptotic behavior of the MLE and MMLE can be studied to justify the usage of the PTP when comparing a misspecified case against a correctly specified one.Referring to the PTP as θ 0 , the true parameter of the Helmholtz model as ϕ, and completely separately considering a correctly specified case where the true model is the delay model with true parameter θ * as was done in ( 40) and ( 41 denotes convergence in distribution [33], [34].Similarly, the MLE has the property that These relationships elucidate that both estimators (asymptotically) belong to the same distribution family, but have different parameters.
If we now let θ * = θ 0 , we observe that the MLE and MMLE have the same mean, but different variances.Correspondingly, choosing to evaluate the delay model at the PTP means that, in this new and correctly specified scenario, the MLE will (in expectation) produce the same parameter estimate that the MMLE produces in the presence of misspecification.As a consequence, the two estimators differ only in their covariance, allowing for a straightforward comparison between the misspecified and correctly specified cases.This justifies our usage of the PTP in the delay model as was done in Figure 4 and Figure 6.
In contrast to this, if one were to evaluate the correctly specified case at a different parameter, say at ϕ as was done in the works [9]- [11], [16], [17], the MMLE and MLE differ both in mean and covariance and a direct comparison would require the computation of a proper distance or divergence between their respective distributions.We can say that ''A reference scenario can be constructed by treating the misspecified model as if it were the true model.Evaluating this new scenario at the PTP isolates the effect of misspecification to the covariance of the MLE and MMLE, making a comparison between the two straightforward.In ultrasound localization, this means that the location parameters and scattering amplitudes in both scenarios will have the same mean values, but different covariance due to the misspecification.Due to the previous discussion point, neither of these estimators will generally coincide with the ground truth.''Choosing to create a reference scenario based on the assumed, misspecified model evaluated at the PTP has additional consequences which we discuss next.

C. BEHAVIOR OF THE CRB
As a final point, we once again draw attention to Figures 4  and 6.In these figures, the MCRB appears to be close to a vertical shift of the CRB, which begs the question of whether the computational cost of the MCRB is warranted.This phenomenon is directly related to the previous point, where it was stated that the reference scenario was constructed by treating the delay model as a new, correctly specified model and evaluating it at the PTP.Doing so eliminates the misspecification and allows the computation of the classical CRB.
As shown in Definitions II.3 and II.4, the MCRB depends explicitly on the PTP θ 0 , which itself depends on the misspecified and true models f X (x|θ) and p X (x).Recalling Theorem 1 and equations (10) and (11), we highlight that the PTP is obtained by projecting the true model onto the misspecified one.This becomes apparent when observing the nonlinear least squares forms in (37) and (38).In this manner, both the choice of true model and its true parameters will affect the PTP θ 0 .
As a result of computing the CRB based on a reference scenario evaluated at PTP, the behavior of the CRB is dictated not only by the assumed model, but also by the true model evaluated at the true parameter.To illustrate this, we provide one final example in Figure 7 in which a third scenario has been considered.The Born model is now taken as the true model, while the delay model remains as the assumed one.The PTP is computed once again using the same ground truth parameters ϕ as for the Helmholtz model in Figure 4c.The Helmholtz and Born models superficially appear to have the same parameter space.The usage of a different model true model, however, introduces a parameter space difference and yields a different PTP.
We observe that the classical CRB for the delay model can behave either as the MCRB using the Helmholtz model or the MCRB using the Born model, depending on which of these is used as the true model from which the PTP is obtained.As such, it is necessary to generate data based on the Helmholtz model to obtain ''realistic'' CRB values in the first place, though they are nevertheless inapplicable due to the misspecification.In particular, if the true model is the Born model, the resulting PTP and classical CRB erroneously convey the message that the estimation task remains simple even when the defects are closely spaced.This is not the case when the true model contains multiple scattering, as shown by the curves labelled ''MCRB H.'' and ''CRB H.'' in Figure 7.
Note that once the data has been generated based on the Helmholtz model and the MMLE has been employed, the computational cost of the CRB and the MCRB is comparable.It is therefore preferable to simply use the MCRB instead.Importantly, experimental design relies on accurately predicting the achievable performance in challenging scenarios, and so the MCRB should be employed over the CRB when misspecification is present.
This can be summarized lightheartedly as ''Constructing reference scenarios to evaluate the impact of model mismatch is a treacherous task in which seemingly innocuous choices regarding the parameters can make a direct comparison of models and bounds difficult or impossible.Parameter differences due to model mismatch and the statistical properties of the MLE and MMLE must be considered.In a poorly constructed reference scenario, neither the expected value nor the covariance of the corresponding MLE task convey any information about the misspecified estimation task.Constructing a reference scenario by evaluating the misspecified model at the PTP is a good choice since, in addition to granting the MLE and MMLE the same expected value, it makes their covariance behave similarly.Nevertheless, in the presence of misspecification, the CRB is not a valid lower bound for any parameter θ of the assumed model.A proper evaluation of the impact of multiple scattering on ultrasound localization tasks should consider all of these factors.''

VII. CONCLUSION
In this work, we have employed the Misspecified CRB to study the impact of the misspecification of multiple scattering in localization tasks where the measurement data stems from a scalar wave phenomenon.We have employed the generalized Slepian formulae to study simulated ultrasound tomography data with the goal of contributing to the discussion on whether multiple scattering is beneficial or detrimental in localization tasks.It was observed that the achievable scatterer localization performance under misspecification of multiple scattering is better than would be predicted by the classical CRB.However, this statement is not to be taken lightly.
The discussion presented in this work aims to highlight the difficulty and nuance involved in a proper evaluation of the impact of multiple scattering in localization tasks.The construction of useful reference scenarios and comparisons is especially daunting, since the influence of parameter space differences and the statistical properties of the MLE and MMLE must be considered simultaneously.To exemplify this, we have recontextualized earlier studies on the impact of multiple scattering through the lens of the MCRB.In doing so, it becomes clear that errors previously attributed to bias are instead caused by parameter space differences between mismatched models.Through this distinction, the MCRB is able to better isolate and describe the influence of misspecification on the localization task than the classical CRB.We advocate for the usage of the PTP in the construction of reference scenarios without model mismatch, as it endows the MLE with desirable statistical properties that account for this parameter space difference.

FIGURE 1 :
FIGURE 1:Illustration of the usage of contrast sources.In Figure1a, an illustration of a generic tomography setup is shown.A transducer ring is placed inside a medium that extends infinitely in all directions.Within the medium, a darker gray color represents a wavenumber that differs from that of the background.The total field within the medium can be expressed in terms of two distinct components.The incident field is shown in Figure1b, where it is clear that the inhomogeneity is no longer present and the excitation produced by a transducer propagates in the homogeneous medium.In Figure1c, the inhomogeneities are shown behaving as sources producing a field that also propagates through the homogeneous background medium characterized by k0.

FIGURE 2 :
FIGURE 2: Scattered fields simulated via the three different models in the presence of two scatterers.Good agreement among the three is seen near the dominant echoes; however, neither the Born approximation nor the delay model exhibit multiple scattering.

FIGURE 4 :
FIGURE 4: Comparison between EMMSE and MCRB for different simulation parameters.The EMMSE is computed for 200 noise realizations when c = 1500 m s −1 and 500 realizations for c = 5000 m s −1 .Figures4a and 4bshow the empirical bias-covariance decomposition of the EMMSE in logarithmic scale.The empirical bias is an order of magnitude smaller than the covariance.Since the EMMSE Q-value is obtained through the Pythagorean addition of the empirical bias and covariance, the influence of the bias is negligible and the covariance coincides with the EMMSE.The EMMSE and covariance additionally coincide with the MCRB.In Figures4c and 4d, the agreement between EMMSE and MCRB is illustrated in linear scale.Even though more noise realizations are used in the computation of Figure4d, the EMMSE appears noisier than in Figure4cdue to the lower SoS contrast and consequently lower SNR.The classical CRB is shown as a reference, and the EMMSE is observed to follow the MCRB instead of the CRB.

FIGURE 5 :
FIGURE 5: Decomposition of the Empirical MSE based on (32).Since the empirical bias is small as illustrated in Figure 4, agreement between the curves means that the error covariance coincides with the MCRB.The small Q-values compared to the true parameter values illustrates that the PTS θ0 of the model ignoring multiple scattering still provides valuable localization information.

FIGURE 6 :
FIGURE 6: Illustration of the Q-images for the MCRB (top row) and CRB (bottom row) for 8 different simulation scenarios.Each image corresponds to a fixed true model (Helmholtz or delay) and a defect speed of sound c with units m s −1 .The assumed model during estimation is the delay model in all cases.The color bar is shown in µm.The values of the MCRB Q-images are lower than those of the classical CRB.The value of the center pixel is undefined, as the two scatterers are at the exact same location.

FIGURE 7 :
FIGURE 7: Illustration of the MCRB and CRB when the true model is the Helmholtz (H.) or Born (B.) model.The classical CRB, where the true model is the delay model, is shown to behave similarly to MCRB B. and MCRB H. simply by changing its parameters to the corresponding PTP.