Echo statistics associated with discrete scatterers: A tutorial on physics-based methods

The software used to produce all predictions in the figures is in the supplementary material at https://doi.org/10.1121/1.5052255. The software is also stored online (Lee and Baik, 2018), where it is subject to future revisions.

Echoes, as measured through the receiver of an active monostatic sensor system, will typically fluctuate from ping to ping as the beam emitted from the system scans across a volume containing scatterers or as the scatterers in that volume move through the beam (Fig. 1). It is essential to understand the echo statistics for accurate interpretation of the scattering and for modeling system performance. Toward that goal, understanding the underlying physical processes that give rise to the fluctuations allows one to accurately predict and interpret the echo statistics. In the simplest case in which the propagation medium is homogeneous and there are no boundaries present (i.e., a “direct path” geometry), the fluctuations are due to a combination of several statistical processes: the interference between overlapping echoes when multiple scatterers are present, the random nature of the echoes from individual scatterers (not including beampattern effects), and the modulation of the echo due to the random location of the scatterer in the beam. Once the geometry is further complicated by heterogeneities in the medium and/or the presence of boundaries, the propagated signals will become refracted and/or rescattered, giving rise to more propagation paths that are potentially random and, in turn, also contributing to the fluctuations. This tutorial presents key concepts and formulations associated with predicting the echo fluctuations over a wide range of scenarios in terms of the physics of the scattering, system parameters, and signal theory.

The statistical behavior of echoes is important across a diverse range of active sensor systems and applications involving the use of either acoustic waves (such as with sonar, medical ultrasonics, or non-destructive testing) or electromagnetic waves (such as with radar or light) to study individual discrete scatterers, assemblages of scatterers, or rough interfaces that cause scattering. Understanding echo statistics has been integral in interpreting radar clutter (Watts and Ward, 2010) and sonar reverberation and clutter (Ol'shevskii, 1978; Gallaudet and de Moustier, 2003; Abraham and Lyons, 2010), sonar classification of marine life and objects on the seafloor (Stanton and Clay, 1986; Medwin and Clay, 1998), medical ultrasound classification of human tissue (Eltoft, 2006; Destrempes and Cloutier, 2010, 2013; Oelze and Mamou, 2016), and non-destructive ultrasound testing of materials (Li et al., 1992). Within each of these areas, echo statistics are used in the detection and classification of scatterers, discriminating between scatterers of interest from clutter (i.e., unwanted echoes), estimating numerical density of scatterers, and determining the performance of sensor systems for use in detection and classification of scatterers. The studies of clutter involve characterizing the statistics of unwanted echoes from features or objects in the environment that have properties resembling the target of interest. The clutter may be due to the sea surface (radar/sonar/laser), marine life or seafloor (sonar), human tissue (medical ultrasound), or grains in solids (non-destructive testing).

The above applications have many elements of statistical theory in common. Those common elements are treated formally in Goodman (1985) and Jakeman and Ridley (2006). There are also notable differences between interpreting echoes from acoustic and electromagnetic systems, such as the presence of shear waves and polarization, respectively, which are summarized in Le Chevalier (2002). Differences in echo statistics between scalar fields (both acoustics and electromagnetic, ignoring shear waves and polarization, respectively) and those fields with polarization effects (electromagnetic only) are summarized in Jakeman and Ridley (2006).

While there has been much work conducted in the area of echo statistics, the focus has generally involved describing the variability of echoes through use of generic statistical functions whose parameters need to be determined from experimental data. Here, “generic” refers to those functions generally devoid of a physical basis and derived solely from random signal theory, such as a Gaussian process. These generic functions include the Rayleigh, Rice, K, Weibull, log normal, and Nakagami-m distributions (App). Only under a narrow range of conditions can some of the generic functions be connected to the physics. There has been a much smaller set of studies in which the statistics were modeled from first principles, and hence being predictive, through a combination of the physics of scattering, sensor system parameters, and signal theory. To date, there has been no text or tutorial with a comprehensive coverage of physics-based phenomena.

This tutorial presents the fundamentals of echo statistics associated with an active monostatic sensor system, written in a manner that is readable by a wide audience. While much of the material is motivated by the use of sonar in underwater applications, the material is broadly applicable to other types of active acoustic systems used in various media including air, solids, and biological tissue, as well as to radar and lasers. Because of the generality of the treatment, the system used to transmit waves and receive the echoes will heretofore be referred to as a “sensor system” regardless of application and whether it involves acoustic or electromagnetic waves. The formulations describe scalar fields, which can be directly applied to either acoustic or electromagnetic signals in the absence of the presence of shear waves or polarization effects, respectively. Only brief reference is made to those latter effects.

The emphasis in this tutorial is on describing the statistics of the magnitude of the complex echo (or echo envelope) in terms of the physics of the scattering, sensor system parameters, and signal processing. More specifically, the emphasis is on physics-based methods, in contrast to the generic approaches referenced above. The methods are focused solely on “first-order statistics” which concern statistical properties of the signal at a single instant in time. Also, the material can apply to single beam (fixed or scanning) and multi-beam systems where the scatterer(s) in each case are randomly located in each beam.

Since there are too many combinations of types of systems, signal processing, and environments to adequately describe in a single paper, the material is focused on the fundamental aspects of echo statistics due to discrete scatterers that are not specific to any particular system, signal processing, or environment. The range of beamwidths, types of scatterers, and number of scatterers in the predictions illustrate the corresponding range of statistical behavior of the echoes from a wide range of system parameters and distributions of scatterers. However, the formulations are completely limited to (1) first-order statistics as described above, and are mostly limited to (2) narrowband signals that are long enough so that the echoes from the scatterers overlap significantly, (3) direct path geometries where there are no boundaries present, and (4) a homogeneous medium. Beyond those limited scenarios, examples are presented of more complex cases involving pulsed signals in which the echoes from the scatterers only partially overlap and the presence of boundaries and/or heterogeneities, including waveguide effects. There are also many more important cases not covered explicitly, but that can be described using these formulations as a basis, including advanced signal processing and complex environments.

A common theme of the material involves predicting the degree to which the statistics of the echo magnitude deviate from the classic Rayleigh distribution (described in more detail later), with a focus on the “tail” of the echo distributions. Experimental observations of echo statistics, particularly in the limit of a large number of scatterers, tend toward the Rayleigh distribution and deviations from that distribution contain information on the scatterers. The tail corresponds to the highest values of echo magnitude and will typically be the part of the echo that is detected above the background noise or reverberation. This tail is shown to contain valuable information as it is a function of the numerical density of scatterers, type of scatterer, bandwidth, and beamwidth.

The tutorial is organized as follows: Secs. II–IV provide qualitative descriptions and illustrations of concepts that are important to echo statistics, a summary of the range of applications that exploit echo statistics, and equations specific to scattering and random processes (not specific to scattering) that are common to many formulations in echo statistics. Sections V–VIII draw from the material in Secs. II–IV in a presentation of echo statistics formulations for a wide range of important physical scenarios. These scenarios include beampattern effects associated with main lobes of various width, narrowband and broadband signals, completely overlapping echoes (long signal) and partially overlapping echoes (short signal), single scatterers and mixed assemblages of scatterers, elongated and randomized scatterers, and geometries involving either a direct path and a homogeneous medium or ones involving the presence of boundaries and/or heterogeneities. Sections VI and VII include a progression of cases, all involving a direct path geometry with a homogeneous medium using narrowband signals that are long enough so that the echoes from the scatterers overlap significantly. The two sections begin with the case involving no beampattern effects (omni-directional beam) (Sec. VI), and then later incorporating beampattern effects (directional beam) (Sec. VII). Section V, besides providing an overview of the material of the rest of the tutorial, outlines how the results in Secs. VI and VII can be extended to more complex cases, such as the ones given in Sec. VIII. Section VIII draws from Secs. VI and VII for selected complex cases involving partially overlapping echoes and the presence of boundaries and/or heterogeneities. Although many details of the derivations are not given in Secs. VI–VIII, references to the concepts and equations in Secs. II–IV are given, as well as references to previously published papers.

Fundamental to the echo statistics are the (1) various types of interference between overlapping echoes from multiple scattering features from a single scatterer, multiple scatterers, and/or multiple propagation paths in a heterogeneous medium or a medium containing boundaries, (2) stochastic nature of scattering by individual scatterers (not including sensor system effects), and (3) random modulation due to random location of the scatterers in the sensor beampattern (Fig. 1). These elements are first discussed qualitatively with important aspects illustrated. Formal mathematical treatments based on the physical processes, sensor system, and signal types are presented later. All processes are described quantitatively in terms of the probability density function (PDF), which is the probability of occurrence of a random variable for each of its values. The PDF and related quantities are defined formally in Sec. IV B 1.

Echoes from one or more scatterers can be decomposed into the sum of the echoes from multiple scattering highlights from each individual scatterer, echoes from multiple scatterers when more than one scatterer is present, and echoes from multiple propagation paths when the medium is heterogeneous and/or there are one or more boundaries present. Regardless of scenario, because of randomness generally associated with the scatterer or environment, the echoes will tend to interfere randomly with each other, causing fluctuations from ping to ping.

In order to understand this interference phenomenon, we first examine the simple case of a single sinusoidal signal whose amplitude (i.e., envelope) does not vary in time. This signal corresponds to the echo from a point scatterer fixed in the beam of a narrowband system. In this case, the amplitude of the signal is constant [Fig. 2(a)]. The PDF describing this single-valued quantity is the delta function [Fig. 2(a)]. Note that when expressing the signal as a complex variable (which is the case for most of this tutorial), then it is the magnitude of the signal that is of interest which is equal to the amplitude of the sine wave in this example.

The statistics change dramatically once another signal with the same magnitude but of random phase is added. This second signal corresponds to the echo from a second identical scatterer (or scattering feature or ray path) whose echo overlaps with that of the first. Now, the PDF of the magnitude of the sum is spread over a range of values, depending on whether the signals added constructively, destructively, or something in between [Fig. 2(b)]. As more signals of equal magnitude and random phase are added, the PDF continues to change and, in the limit of an infinite number of signals of random phase, the PDF converges to the Rayleigh PDF [Fig. 2(c)].

This example involved single-frequency signals of infinite extent. This idealized signal can be used to approximate pulsed (or “gated”) sine wave signals used in many sensor systems where echoes from scatterers significantly overlap. The same principles also apply, although with more complexity, to pulsed signals (narrowband and broadband) of much shorter extent where the multiple signals generally overlap only partially as presented in Sec. VIII A.

The simplest of scatterers is a point scatterer, as it gives the same echo regardless of orientation. In this case in which effects from the sensor system are ignored (i.e., omnidirectional beam) and when using a sinusoidal signal, the echo magnitude from a single point scatterer is constant with a delta function PDF [Fig. 3(a)]. For more realistic scatterers of finite extent, such as elongated ones whose orientation and/or shape may change in time, the echo magnitude will vary from ping to ping. Now, the PDF is significantly broadened and the echo statistics are more complex [Fig. 3(b)].

There is some correspondence between this and the former case involving sinusoidal signals. The echo from a point scatterer corresponds to a single sine wave with a constant magnitude. In some cases, an echo from an irregular finite-sized scatterer can be decomposed into echoes from different parts of the body of the scatterer, which corresponds to the addition of multiple random-phased sinusoidal signals. Also, there are conditions under which the echo PDF from a complex scatterer can be the Rayleigh PDF, which corresponds to the sum of an infinite number of randomly phased sine waves.

The echo measured through the receiver of the sensor system depends, in part, on the location of the scatterer in the two-way beampattern. The beampattern will modulate the echo according to the location, with stronger values being associated with the scatterer being near the center of the main lobe and weaker values corresponding to locations well away from the center. For a randomly located scatterer, the modulation will correspondingly be random, adding to the variability of the echo. The distribution associated with this random modulation is referred to as the “beampattern PDF.”

The randomizing effects of the beampattern on the echo from a randomly located object can be illustrated by varying the solid angle within which the object is allowed to move (Fig. 4). In this simple example, a circular piston transducer is used to both send and receive the signal and the solid angle is centered about the axis of symmetry of the transducer (which is the center of the main lobe). A point scatterer is allowed to move randomly within the solid angle at a fixed distance from the transducer. In the first case, the solid angle is 0 and the scatterer is restricted to remain fixed in the center of the beam where the response is maximum [Fig. 4(a)]. The corresponding echo is single valued from ping to ping and its PDF is the delta function [right panel of Fig. 4(a)]. As the solid angle increases from zero allowing the scatterer to randomly move a greater amount across the beam, the ping-to-ping variability in the echo increases correspondingly [Figs. 4(b)–4(d)]. When only the portion of the main lobe above the highest sidelobe is involved in the motion [Figs. 4(a)–4(c)], then the echo PDF is either single valued or approximately power-law-distributed. However, once sidelobes are involved, then the echo PDF is more complex, as a non-monotonic characteristic is present due to the sidelobe structure [“spiky” section of curve in Fig. 4(d), where values of echo are lower].

These effects will be shown in Sec. VII A 4 to be qualitatively similar over a wide range of beamwidths. In essence, regardless of beamwidth (several degrees or several tens of degrees), if the scatterer is randomly located within a solid angle containing at least most of the main lobe of the beampattern, then the echo will correspondingly be modulated over a wide range of values of the beampattern.

Understanding and quantifying the statistical variability of a signal is useful in a diverse range of applications. Below are described two common uses of echo statistics, one in which the variability of the observed signal is used to infer important information regarding the scatterers, and the other in which the degree to which the signal varies for a given scenario is predicted to understand the error or uncertainty regarding the expected signal. Specific examples are given in Sec. III A where echo statistics are used as an inference tool spanning use of sonar, radar, and medical ultrasound.

Regardless of the field and type of sensor system, one can exploit properties of the variability of the echo to make inferences of important quantities such as number of scatterers, scatterer characteristics, discriminating between and classifying different types of scatterers, and probability of false alarm when detecting a scatterer of interest that is interspersed with other unwanted scatterers. Key to the success of the inference is using a model of the echo statistics whose parameters can be linked to the physical properties of interest. The description of such models is the focus of this tutorial.

Once the statistics model is determined for the particular application, there are numerous methods to infer parameters of modeled echo statistics from data including use of least squares, maximum likelihood estimators (MLE) (Azzalini, 1996), and method of moment estimators (MME) (Joughin et al., 1993). Performance of an unbiased estimator can be evaluated formally through the Cramer-Rao lower bound (CRLB) method which estimates a lower bound of the variance of the estimator (Hogg and Craig, 1978). Abraham and Lyons (2002b) provides an example of use of MLE, MME, and CRLB for inference of modeled parameters from data (Appendix B of that paper).

Many applications using echo statistics as an inference tool are described in citations given in Sec. I. Selected examples spanning three types of sensor systems—sonar, radar, and medical ultrasound—are briefly summarized below.

The shape of the echo PDF can be used to infer the number of scatterers. This topic has especially seen much attention spanning many applications, with attempts to connect parameters of either physics-based or generic PDFs, such as the K PDF, to numerical density of scatterers. Stanton et al. (2015) used physics-based methods (such as described in this tutorial in Sec. VII B) to relate the shape of the PDF of the magnitude of the echo to the number of scatterers in a sonar beam in the ocean. The inferences were conducted for several different types of scatterers (“bottom-like,” “compact stationary,” and “compact non-stationary”) that ranged in densities from resolved to unresolved. Lee and Stanton (2015) used similar physics-based methods with a sonar beam to infer the numerical density of fish (more details on this broadband method given in Sec. VIII A). Other studies have related a parameter of a generic PDF to numerical density of scatterers. For example, Abraham and Lyons (2002b) analytically related the shape parameter of the K PDF to numerical density of scatterers (for certain conditions; general and not specific to sonar) and applied the results to using sonar to estimate numerical density of scatterers on the seafloor. Tunis et al. (2005) performed a controlled laboratory experiment with various dilute solutions containing cancerous cells (acute myeloid leukemia and prostate adenocarcinoma) to empirically relate a parameter of the gamma distribution of the medical ultrasound echo to numerical density of the cells.

For the general case in which there is a mixture of different sized scatterers, Lee and Stanton (2014) have formulated the echo PDF based on various physical parameters—sensor beampattern, scattering amplitudes, and numerical density of each type of scatterer. This method is not specific to any type of sensor system (and is described in Sec. VII C) and the simulations demonstrated conditions under which the number of scatterers (especially within the type that dominated the echo) could be inferred (see Figs. 6 and 7 of that paper).

For a single scatterer of interest interspersed within an aggregation of other unwanted scatterers (in a volume or on a surface), echoes from the scatterer of interest can be discriminated from those containing only the surrounding unwanted or “background” scatterers. Here, the echo containing the scatterer of interest will also be contaminated with echoes from the background. Ferrara et al. (2011) conducted an experimental study on the ocean involving use of radar to detect and classify echoes from ships and oil rigs. In their data, there were echoes from two types of regions containing (1) both the scatterers of interest (ships or oil rigs; one at a time) and the sea surface and (2) only the sea surface. The experimental radar-echo statistics data were fit to the generalized K PDF. They demonstrated that combinations of generalized K parameters (such as ratios of the parameters) could be used to unambiguously distinguish between echoes involving the scatterers of interest (ships and oil rigs) and the background (sea surface) scattering from those due to the surrounding sea surface alone. There was only one false alarm out of 229 detected targets.

While the previous case in Sec. III A 2 involves classifying a scatterer of interest when its echo is confounded with that of the surrounding background, the case simplifies once the echo is resolved from the background. This can occur when the scatterer is suspended in the volume away from any boundary or other neighboring scatterer. In this case, through inversion methods, variability due to the scatterer's random location in the beampattern can be removed to isolate echo variability due to the target alone. The target echo variability (with no beampattern effects) contains information on the scatterer, such as its size and degree to which it is elongated or rough.

For example, in Clay (1983), the PDF of the echo (as described in Sec. VII A 1 of this tutorial), which is a function of both the beampattern PDF and the PDF of the scattering amplitude of the scatterer, is formulated in terms of a convolution involving each of those individual PDFs. Through knowledge of the beampattern properties, the effects of the beampattern are removed from the echoes through deconvolution, leaving only the PDF of the scattering amplitude. The statistics of the scattering amplitude are then used to classify the scatterer in terms of its size and type. Although the paper was intended for use of sonar to detect and classify echoes from fish, the equations are general and applicable to any sensor system. Clay's method was refined in Stanton and Clay (1986).

Oelze and Mamou (2016) reviewed a number of previously published studies concerning use of medical ultrasound to detect and classify in vivo soft tissue. One such study (Fig. 3 of that paper) combined the spectral content of the broadband echo [to estimate effective scatterer diameter (ESD)] with echo statistics data fit to the homodyned K distribution. The feature analysis plot of ESD versus two parameters of the homodyned K distribution provided strong discrimination (with no overlap) between echoes from three different types of cancerous tumors (mammary) in rodents (rat fibroadenomas, mouse carcinomas, and mouse sarcomas).

In order for the echo statistics to best be exploited, there are further considerations in using echo statistics as a tool with a practical sensor system.

The degree to which the echoes are non-Rayleigh contains valuable information such as numbers of scatterers that can be inferred, as summarized above. Once the echoes are Rayleigh distributed, the amount of information to be inferred is limited. The beampattern causes echoes that may be otherwise Rayleigh-like to be non-Rayleigh through reducing the effective number of scatterers that are “seen” by the narrow mainlobe of the system. Furthermore, echoes that are non-Rayleigh before beampattern effects will deviate even more from the Rayleigh PDF once the beampattern is accounted for. The process of inferring information from the field of scatterers can therefore be enhanced or even optimized through varying, when possible, the width of the main lobe of the beampattern so that the echoes will be non-Rayleigh over the range of expected values.

In any practical sensor system, there will be various sources of noise (including electronic noise) and background reverberation that will normally dominate the lower values of echoes. Thus, any information to be inferred must generally come from the higher portion of the echo PDF—that is, the “tail.” Many analyses focus solely on properties of the tail that is above the detection threshold of the system.

In many applications, it is desired to predict the magnitude of the signal after it has propagated through an environment. The applications can deal with ones such as those described above in which one wishes to infer properties of the environment or scatterers from the signal or other types of applications in which the performance of the sensor system is being studied. Ideally, one wants to know all key parameters of the environment and scatterers so that the properties of the signal that propagates through the environment can be known to 100% certainty. However, due to a combination of the uncertainties of the values of the many parameters of the environment [e.g., roughness of surface(s), material properties, etc.] and scatterers (e.g., number, type, orientation, etc.) as well as random variations of those same quantities due to naturally occurring processes, the predicted signal in a natural environment is normally a random variable.

To quantify the degree to which the signal varies (that is, predict what is commonly called its error or uncertainty), the random variables associated with parameters of the environment and scatterers must first be estimated. Using these random variables, the variability of the signal can then be predicted through use of physics-based equations as described below. Regardless of the approach used, there are different degrees of variability. First, predicting the variability of the signal due to random variations in the environment and scatterers from the expected mean properties. A higher order variability would be due to errors in those estimates of the random variations in the environment and scatterers and associated deviations from the predicted variability of the signal (that is, an error in a prediction of a signal PDF).

Beyond the above considerations in which there is no noise in the system, there is also uncertainty in the predictions associated with random noise. This noise can be due to a combination of ambient noise in the environment and electrical system noise. Also, diffuse background reverberation is sometimes considered to be part of the noise. Although this topic is outside the scope of this tutorial, it is an important consideration in modeling practical systems. The noise can be accounted for through various methods which include (1) modeling the noise as an additive random complex variable to the (complex) echo signal (Jones et al., 2017) and (2) modeling the echo magnitude PDF as a “mixture” PDF, which is the sum of the noise-free echo magnitude PDF and the noise-only PDF (Stanton and Chu, 2010; Abraham et al., 2011).

Fundamental equations are first given describing deterministic scattering such as for a fixed location or orientation of the scatterer. In order to describe the scattering for a randomized case such as when the location or orientation vary randomly, general equations describing statistical processes that occur in such scenarios are given, which are then incorporated into the equations for deterministic scattering to describe randomized scattering. All of the below equations relate the echo and its fluctuations directly to the physical properties of the scatterers and sensor system. This type of treatment is referred to as “physics based.”

The scattering described in this section is from a single ping or statistical realization. This deterministic description will then be randomized for use in predicting the statistics of the echoes. The scattering geometry involves use of an active monostatic sensor system where the transmitter and receiver are collocated (Fig. 1). The scattering can involve one or multiple scatterers. Also, all analyses involve scatterers distributed within a thin shell of constant radius with the sensor system at the center. The thickness of the shell is much smaller than its radius.

This (shell) geometry is for modeling a long narrowband signal in which echoes from all of the scatterers within the shell completely overlap. For a sensor system that uses short pulses to detect a more broadly distributed set of scatterers (such as in many practical applications), the thickness of this thin shell can be related approximately to the duration of the pulse. A rigorous approach to modeling pulsed systems where the echoes generally only partially overlap is given in Sec. VIII A.

The voltage V_s received by the sensor system due to the echo from a single scatterer is

where $VT$ is the voltage applied to the transducer, $GT$ is the transmitter response (conversion factor of applied voltage to acoustic or electromagnetic field) which is equal to the acoustic or electromagnetic signal, respectively, at reference distance $r0$ per unit applied voltage $VT$⁠, $GR$ is the receiver sensitivity (conversion of acoustic or electromagnetic field to voltage) which is equal to the voltage signal at the output of the transducer per unit acoustic or electromagnetic signal, respectively, incident at the transducer, $r$ is the distance between the transducer and scatterer, j = $−1$⁠, $ω$ is the angular frequency of the sinusoidal signal, $k$ is the acoustic or electromagnetic wavenumber ( = 2$π$/λ, where λ is the wavelength), α is the absorption coefficient of the medium so that $e−2αr$ is the two-way loss due to absorption, and $b(θ,ϕ)$ is the two-way beampattern of the sensor system whose values lie in the range [0,1]. The term $b(θ,ϕ)$ is the product of the beampatterns of the transmitter and receiver and the terms $θ$ and $ϕ$ are the angular coordinates of the scatterer. Specifically, $b(θ,ϕ)=bT(θ,ϕ)br(θ,ϕ)$⁠, where $bT$ and $br$ are the transmitter and receiver beampatterns, respectively. The term $fbs$ is the backscattering amplitude of the scatterer and is a complex variable.

In this formulation, the signal is assumed to be at a single frequency (i.e., narrowband) of infinite temporal extent. Also, the acoustic and electromagnetic fields associated with $GT$ and $GR$ are assumed to be scalar quantities. For acoustics, this scalar quantity is pressure, the compressional component of the field, and assumes no shear component in the medium (although conversion of compression to shear wave can take place within the scatterer). For the electromagnetic field, the scalar quantity is one polarization component of the electric or magnetic field. This assumption in this latter case treats each polarization component independently and ignores coupling between the components [Chap. 4 of Goodman (1985); Chap. 4 of Jakeman and Ridley (2006)].

The target strength of the scatterer can be expressed in terms of the backscattering amplitude and differential backscattering cross section $σbs$ as

where $σbs$ = $|fbs|2$ and the units of $fbs$ (m) and $σbs$ (m²) are suppressed. Note that the term $σbs$ should not be confused with a similar representation for backscattering cross section, $σ$⁠, that is commonly used where $σ=4πσbs$ and TS = $10 log(σ /4π).$

For simplicity in the analysis, all parameters of the system and measurement are assumed to be constant and will be suppressed in the following equation. From Eq. (1), the magnitude of the echo voltage due to a single scatterer as received through the monostatic active sensor system is now given by

where $VT$⁠, $GT$⁠, $GR$⁠, $r,$ and α in Eq. (1) have been suppressed and $r0=1$m. Here, the magnitude of the signal is based simply on the absolute value of the signal.

The voltage received by the sensor system due to the echo from an aggregation of N scatterers is

where $ri$⁠, $fbs(i)$⁠, and $(θi, ϕi)$ are the range, backscattering amplitude, and angular location of the ith scatterer, respectively. As with Eq. (1), the signal is assumed to be at a single frequency of infinite temporal extent. With this type of signal, there is 100% overlap between the echoes from all individual scatterers. The simple summation of echoes from individuals reflects the assumption that only single-order scattering is being considered and higher-order scattering (i.e., re-scattering of echoes between individuals) is assumed to be negligible.

From Eq. (5), and suppressing the constants of the system and measurement in a manner similar to that with the individual scatterer described above, the magnitude of the echo voltage due to N scatterers as received through the sensor system is given by

where the magnitude of the echo voltage from the ith scatterer as received through the sensor system is

The assumption that all parameters of the sensor system and measurement are constant requires the range from the sensor system to each scatterer to be approximately the same so that the differences in losses due to spreading and absorption in the $ri−2$ and $e−2αri$ terms, respectively, in Eq. (5) associated with each scatterer are negligible (i.e., $ri−2∼r−2$ = constant and $e−2αri$ ∼ $e−2αr$ = constant). Specifically, all scatterers are assumed to be located within a thin shell whose thickness is small compared with the radius of the shell. However, minor differences in range within the shell can still lead to significant changes in phase of the echo from individual scatterers, especially when the acoustic or electromagnetic wavelength is comparable to or smaller than the shell thickness. For simplicity in the formulation, all phase shifts associated with the ith scatterer are in the term $Δi$⁠, which includes phase shifts due to differences in range within the shell (i.e., 2kr_i) and due to the scatterer and beamformer. For acoustic/electromagnetic wavelengths that are small compared with the shell thickness, $Δi$ will generally vary randomly and uniformly in the range [0 2$π$] for randomly distributed scatterers.

The statistics described herein involve the probability of occurrence of random variables. This is in contrast to other types of statistics such as statistical tests (e.g., the t- and Kolmogorov-Smirnov tests). The principal statistical quantities used are the probability density function (PDF or p), cumulative distribution function (CDF), and probability of false alarm (PFA), which are interrelated from the following expressions for one-dimensional continuous random variables (Ol'shevskii, 1978; Papoulis, 1991; and Goodman, 1985). Expressions involving multi-dimensional random variables will appear later in context. While the below expressions are general, the integration limits associated with echo magnitude statistics later will reflect the fact that the magnitude (x) is always positive and $p=0$ for $x<0$⁠.

Note that the CDF is rigorously referred to as the “distribution” or “distribution function.” However, the PDF is also commonly referred to as a “distribution” as well as “frequency function” throughout the literature. While “PDF,” “distribution,” and “distribution function” are commonly interchanged with no change in meaning (in context), the term “PDF” will be principally used herein.

The infinitesimal probability dP_X of a random variable X occurring in the differential interval [x, x + dx] is expressed in terms of the probability density function $pX(x)$ of X,

For a finite interval $[a,b]$⁠, the probability is calculated through the following integral:

Once $a$ and $b$ are extended to $−∞$ and $+∞$⁠, respectively, the integral over all values of X is equal to unity.

The probability of a random variable occurring for all values of X up to an arbitrary point x is determined from the integral in Eq. (9) and is referred to as the cumulative distribution function,

As stated above, once these equations are applied to echo magnitude statistics, $pX(u)=0$ for $u<0$⁠, thus the lower limit of this integral would be zero.

The PDF can be determined from the CDF simply by taking the derivative of the above expression,

The probability of false alarm PFA is commonly used to determine the probability of occurrence of a random variable occurring for any value higher than an arbitrary value (Chap. 7 of Ol'shevskii, 1978). From Eq. (9), the PFA of X for all values above x can be related to the integral of the probability density function as

Using the fact that the integral of the PDF over all values of its argument is equal to unity, Eqs. (10) and (12) can be used to express the PFA in terms of the CDF,

The PFA presented in this tutorial is mathematically equivalent to the probability of detection (PD). The nomenclature varies depending on the context of application. For example, when the PDF is used to describe the unwanted background echoes or “noise,” the PFA gives a measure of the probability that the source of scattering is not from the target of interest. Conversely, when the PDF is used to describe the anticipated echo from a scatterer (or “target”) of interest, the PD gives a measure of the probability of detecting the target of interest. The limit, x, in Eq. (12) is the threshold above which the PFA and PD are calculated (Chap. 7 of Ol'shevskii, 1978).

Some PDFs are, conveniently, closed-form analytical solutions that can be calculated directly. For example, the Rayleigh and K PDFs given later, as well as the Gaussian PDF are in a simple analytical form and are straight forward to calculate. Other analytical solutions, such as those given in integral form, are also straight forward to calculate through numerical integration. However, PDFs for many realistic cases are generally not in closed form and require numerical simulations involving the scattering equations to estimate them. For example, once parameters of the sensor system and scatterers are specified, models of the echo at the signal level (not PDF level) such as in Eq. (6) are simulated to create the statistics of the echo.

A common method to estimate PDFs is through use of Monte Carlo simulations which generally involves making calculations of many statistical realizations of the signal to form an ensemble from which the PDFs are estimated. For example, in Eqs. (6) and (7), each ith scatterer will have a distribution of scattering amplitudes $fbs(i)$ and locations (⁠$θi, ϕi$⁠). Each realization involves randomly selecting one of those scattering amplitudes and locations for each scatterer, then calculating the signal using Eq. (6). The process is repeated many times, with each selection of random values being statistically independent (when independence is required) from the other selections.

Use of Monte Carlo simulations tends to be more general than analytical methods as there are fewer assumptions and, hence, limitations. For example, the scatterers can be arbitrarily correlated in space (such as being in a thin layer or small patch) and can have arbitrary phase-shift distributions. While many analytical models are restricted to cases involving independently distributed scatterers whose phases are uniformly distributed over [0 2$π$], use of the Monte Carlo simulations do not have such restrictions.

Once the many simulations have been completed, the PDF is commonly estimated by putting the realizations of the signal into “bins” to form a histogram. For an analysis of echo magnitude statistics in the example described above, the result of each calculation is put into a magnitude bin (i.e., quantized value of echo magnitude) so that a histogram can be formed. These simulations require many realizations so that there can be correspondingly many narrow bins in order to produce a histogram that is an accurate representation of the actual PDF.

This binning approach is intuitive and is a method used in this paper, when appropriate, due to its simplicity. Conditions under which this method are used depend upon a combination of the type of structure in the PDF and corresponding number of realizations required to model that structure. In some cases, such as for a smoothly varying PDF, the computation time is reasonable. Calculations of PDFs for other applications where there is structure such as the presence of narrow peaks or nulls in the PDF curve may require more realizations and correspondingly significant computer time. When making too few calculations in this latter case, there can be artifacts in the result, such as smoothed or completely missed peaks or nulls. Thus, when there is the presence of narrow peaks or nulls in the PDF, a closed-form analytical method is used in this paper to determine the PDF, when possible. Beyond these approaches, the kernel density estimation (KDE) method was used to reduce the number of realizations needed to produce a reliable estimate of the echo PDF (Botev et al., 2010; Lee and Stanton, 2015; Scott, 1992). The calculations illustrated in this paper typically involve 10⁷ realizations.

Finally, for applications that extract information from the tail of the PDF, estimation methods such as importance sampling can be used to reduce the variance in the estimate and to increase the efficiency of the Monte Carlo process through selectively sampling the more desired (tail) samples (Agapiou et al., 2017).

Depending upon the types of features one is investigating in a PDF, the curves will either be plotted on a linear-linear or logarithmic-logarithmic scale. While the former scale may be more intuitive, the latter is especially useful when examining the tail of the PDF which typically has low values relative to the maximum. Choice of type of scale influences how the bins are determined. For linear-linear plots of PDFs, equally-spaced bins for the horizontal axis are normally used. However, when plotting PDFs on a logarithmic-logarithmic scale, the width of the bins should be equal on a log scale, which is non-uniform on a linear scale. Otherwise, if the bins were equal on a linear scale, but plotted on a log scale, the density of points on the plots would increase throughout the plot, and not fully characterize the shape of the PDF.

The probability of a variable occurring over any of the values of the random value x over the entire range is, by definition, unity. Therefore, the integral over x of any PDF over all values is unity. PDFs are commonly derived with a constant factor introduced that is determined through normalizing the area under the PDF curve to unity. From this property, it follows that the CDF will begin at a value of 0 for the smallest value of x and reach its maximum value of unity at the largest value of x. Similarly, the PFA will begin at unity and decrease to the value of 0 for the corresponding smallest and largest values of x, respectively.

In some applications, it is also important to normalize the (horizontal) scale associated with the random variable. This can be the case when the calibration of the system is not known accurately, the propagation loss of the signal in the medium is not known accurately, or when only the shape of the PDF, CDF, and PFA are of interest regardless of the echo strength. Through normalization, only the relative values of the random variable will be considered. One convenient approach is to normalize the random variable by its root-mean-square (rms) value $⟨x2⟩1/2$ and plot the PDF, CDF, and PFA versus the random variable divided by $⟨x2⟩1/2$⁠, where $⟨⋯⟩$ is the average over a statistical ensemble of values. In this case, the area under the PDF curve (with an argument normalized by $⟨x2⟩1/2$⁠) is preserved under the transformation and is unity.

Regardless of the type of scale (linear-linear or logarithmic-logarithmic) or uniformity of spacing of bins, all normalizations are first calculated on a linear-linear scale. For example, an equation such as Eq. (9) is on a linear-linear scale and can be used to normalize the PDF to unity while accounting for non-uniform spacing in the integral.

Random fluctuations of echoes involve several fundamental statistical processes. For example, in Eq. (1), the beampattern is shown to be a function of the angular coordinates. In general, the scatterer will be randomly located in the beam, making the angular coordinates of the scatterer random variables. Since the beampattern is a function of the random angular coordinates, then the beampattern function is, in turn, a random variable for a randomly located scatterer. The scattering amplitude in Eq. (1) is also generally a randomly variable due to the random nature of the scatterer. Since the echo $ẽ$ in Eq. (4) as measured through the receiver of the sensor system is the product of the two random variables, the beampattern function and scattering amplitude, then $ẽ$ is also a random variable. Finally, once there are multiple scatterers in the beam of the sensor system, the resultant echo $ẽ$ [Eq. (6)], as measured through the receiver, will be the sum of the random individual echoes [Eq. (7)] and will, in turn, be a random variable.

These statistical processes—function of a random variable(s), multiplication of two random variables, and addition of random variables—are of wide applicability, are not specific to sensor systems or scattering, and appear in standard textbooks on statistics. Formulas summarizing these general processes are given below for later reference in the scattering problem. While only the simplest of cases involving one or two random variables are given, formulas involving more random variables are given in the references and/or later in context of the application.

If the function Z is a function of the random variable X, then Z(X) is also a random variable. The formulations relating the PDF of Z to the PDF of X are based on the fundamental principle that the probability of occurrence of an event in one space (X in this case) is the same as that in the transformed space (Z in this case). The resultant PDF $pZ$(z) is then given by one of two equations depending upon whether Z(X) varies monotonically or non-monotonically with respect to X. Specifically, Z(X) is monotonic with X if it either solely increases or solely decreases over the range of X such that for any value of Z, there is only one (unique) value of X. Conversely, for the non-monotonic case, Z(X) both increases and decreases over the range of X so that there can be multiple values of X for a given value of Z [both of these cases and the below equations are described on pp. 23–27 of Goodman (1985)].

For the case in which Z(X) varies monotonically with respect to X over the entire range of X, then the following expressions can be written where the differential probabilities in the two spaces are equated to each other,

From Eq. (8), this can be expressed in terms of the PDFs of X and Z,

Rearranging terms yields an expression for the PDF of Z in terms of the PDF of X for this monotonic case,

where the absolute value sign is used to keep the expression for the PDF positive.

In the more complex case in which Z(X) varies non-monotonically with respect to X over the range of X, the PDF is described by a similar equation, but summed over M contiguous segments where Z(X) varies monotonically within each segment,

Here, x(z) and x_m(z) in Eqs. (16) and (17) are the inverse functions z⁻¹(x) and z⁻¹(x_m), respectively. In practice, these inverse functions can be determined numerically from the forward analytical function, plots, or tables of z(x) and z(x_m).

The above analysis involving a function of one random variable is extended to the case of a function of two random variables. In this case, if the function Z is a function of the random variables X and Y, then Z(X, Y) is also a random variable. Relating the PDF of Z to the PDF(s) of X and Y involves the same process as in the previous case of one random variable in which the probability of occurrence of an event in one space is set equal to that of the other space. This process generally involves first determining the Jacobian of the transformation relating the two spaces, although that will not be shown explicitly below (Papoulis, 1991).

From Eq. (6-35) of Papoulis (1991), the probability of Z occurring for any value below z is given in terms of x and y as

where P_Z (⁠$Z≤z$⁠) is also the CDF_Z. Here, $pX,Y(x,y)$ is the joint probability density function of the random variables X and Y, and D_Z is the region or regions in the xy plane containing values of x and y where Z(X,Y) $≤z$ (D_Z is illustrated in Fig. 6-7 of Papoulis, 1991). This equation for P_Z (⁠$Z≤z$⁠) is a two-dimensional form of Eq. (10). The PDF of z can be expressed by taking the differential of $PZ$ (⁠$Z≤z$⁠) above,

where dD_z is now the differential region or regions(s) in the xy plane whereby the values of x and y are bounded by the differential area determined by the range $z≤Z≤z+dz$ [Eq. (6-36) and Fig. 6-7 of Papoulis, 1991].

This equation is complex to solve and depends upon the characteristics and form of $pX,Y(x,y).$ For simple forms such as Z = XY and Z = X + Y, where X and Y are independent random variables, then the solution for p_Z(z) in each case is in closed form. Since those two cases are used throughout this tutorial, they are treated separately in Secs. IV C 4 and IV C 5. The more complex case in which Z is a general function of X and Y is used only once later in the tutorial and the solution will be given in the context of that application (Sec. VII A 7).

If X and Y are random variables, then the product Z = XY is also a random variable as described in Sec. IV C 3. If X and Y are independent of each other, then $pX,Y(x,y)$ = p_X(x)p_Y(y) in the integrand in Eq. (19), where p_X(x) and p_Y(y) are the PDFs of the random variables X and Y, respectively. Inserting this product of the two PDFs into the integrand in Eq. (19), the PDF $pZ(z)$ of the product Z = XY can be shown to be

This equation is from Eq. (6-74) in Papoulis (1991). In that book, the equation is derived through a method involving use of a Jacobian of the transformation to map one coordinate system to another. However, this equation can also be derived directly from Eq. (19) of this tutorial [which is Eq. (6-36) of Papoulis (1991)] using the change of variables method illustrated in Papoulis for the ratio of two random variables (p. 138 of the book). Using that method in this case for the product of two random variables (Z = XY), the change in variables y = z/x is used, and the area dxdy for dD_z is mapped to the area (1/|x|)dxdz. Through this mapping, the double integral for dD_z is replaced with a single integral over x. Replacing dxdy in Eq. (19) with (1/|x|)dxdz, the dz drops out of both sides of the equation and the integral is only over dx as shown. The absolute value sign is used for the variable x so that the differential area will be positive for all values of x. Note also that the term |x| in the factor (1/|x|) for the area is equal to the absolute value of the Jacobian of the transformation in the derivation of Eq. (6-74) in Papoulis (1991).

Once Eq. (20) is used in physical applications, the range over which one or more of the physical parameters may be constrained and its corresponding PDF will be zero outside of that range. The integration limit(s) may reflect that constraint by only spanning the range over which the PDF is non-zero as shown later.

There is a variety of methods to evaluate the PDF of the sum of independent random variables, ranging from purely analytical to purely numerical. Sometimes, a “purely” analytical method still requires numerical evaluation, such as when an integral or series summation are involved and numerical integration or summation are performed, respectively. Two common methods are discussed below: the method of characteristic functions and Monte Carlo simulations.

A commonly used analytical method involves use of characteristic functions (CFs) where the CF of a random variable is the Fourier transform of its PDF (Goodman, 1985). Addition of an arbitrary number of independent random variables involves first taking the product of their corresponding CFs. This product is the CF of the sum of the random variables. The PDF of that sum is then the inverse Fourier transform of the CF product.

This CF approach can be derived from Eq. (18) for the case of Z = X + Y where Z, X, and Y are all random variables. Since the random variables, X and Y, are independent of each other, then $pX,Y(x,y)$ = p_X(x)p_Y(y). Using this relationship in the integrand of Eq. (18), the PDF of Z, $pZ(z)$⁠, can be shown to be the convolution of the two functions, p_X(x) and p_Y(y). This convolution can then be shown to be equivalent to the product of the Fourier transforms of p_X(x) and p_Y(y). Since these Fourier transforms are, as defined above, the CFs of the two functions, then the method of characteristic functions follows as described above (Goodman, 1985; Papoulis, 1991). The method is extendable to the sum of an arbitrary number (N) of independent random variables by first expressing the convolution integral by formulating the sum of two random variables where one of the random variables is the sum of N−1 random variables and the other random variable is the remaining variable. The PDF of the sum of the N−1 random variables is determined through a similar process involving the sum of N−2 random variables, and so on. After completing this iterative process, the PDF of the summed N random variables is related to the product of the Fourier transforms of the PDFs of the N random variables.

Acoustic and electromagnetic signals are complex and normally constructed of a real and imaginary term, making them two dimensional. Since the method of characteristic functions is extendable to multi-dimensional variables, this method can be applied to determine the PDF of the sum of complex signals. For the case in which the phase of the summed signal is uniformly distributed [0 2$π$] and each component of the signal has a zero mean (i.e., a “circularly symmetric signal” in the complex plane), then the CF and PDF of the signal magnitude are a Hankel transform pair. Application of the CF to calculate the PDF of the magnitude and magnitude squared of complex signals is summarized in Jakeman and Ridley (2006, Chap. 4), including an extension to signals where the phase is not uniformly distributed. Methods to numerically evaluate the PDF for circularly symmetric signals (via the CF and Hankel transform) are given in Drumheller (1999).

Barakat used a broadly similar approach to the Hankel transform method by extending the 1D CF method to circularly symmetric complex signals through constructing an orthogonal component of the sum, equal in magnitude to the original single component of the sum, resulting in an exact, analytical expression for the PDF of the magnitude of the sum of complex random variables (Barakat, 1974). However, our experience in applying the Barakat approach has resulted in convergence issues due to truncation of the infinite series summation that must be evaluated (Chu and Stanton, 2010; Lee and Stanton, 2014).

The method of Monte Carlo simulations is discussed in more general terms in Sec. IV B 2 and will only be briefly summarized here in the context of adding random variables. This is a commonly used numerical approach that involves simulating a statistical ensemble of a large number of realizations of the process of interest so that the PDF can be formed. Performing Monte Carlo simulations to evaluate the sum of random variables in predictions of scattering will provide a stable and accurate solution (see, for example, Stanton et al., 2015). When using this method to simulate signals associated with a sensor system, random variables are normally added through (1) a phasor addition in the frequency domain of single frequency or narrowband signals of long enough extent that the echoes are completely overlapping [such as with Eq. (6)] or (2) addition in the time domain of short signals when the echoes are only partially overlapping and/or when broadband signals of any duration are used (Sec. VIII A). In the case of phasor addition, the signals are first represented in complex form and then the real and imaginary components are added separately before being recombined to calculate the signal magnitude.

In the limit of the sum of an infinite number of independent complex random variables, drawn from identical distributions with uniformly distributed phases, each of the two independent components of the sum tends to a Gaussian PDF, with zero mean and equal variance. This is referred to as the central limit theorem (CLT) and is integral to many treatments of random variables (Goodman, 1985; Jakeman and Ridley, 2006). The statistics of the magnitude of the sum can be shown to be the Rayleigh PDF,

where $λR=⟨x2⟩$ is the mean square magnitude.

From Eqs. (10) and (13), the CDF and PFA associated with the Rayleigh PDF are

where the lower bound in the integral in Eq. (10) is zero since $pRay(x)=0$ for $x<0$⁠.

The Rayleigh PDF is widely used in describing echo statistics. It is commonly used as the “starting point” in describing the statistics of the echo magnitude, especially when there are many scatterers or many highlights from an individual scatterer contributing to the echo. When the statistics do not follow the Rayleigh PDF, deviations of the statistics from the Rayleigh PDF are frequently described. The deviations in the higher values of the echo magnitude, i.e., the “tail” of the distributions, are of particular interest. The term “non-Rayleigh statistics” is commonly associated with those distributions that deviate from the Rayleigh PDF. The Rayleigh PDF and associated CDF and PFA are illustrated (Fig. 5).

As discussed in Sec. I, various important aspects of echo statistics will now be examined in detail. The treatment will draw from the concepts and equations given in Secs. II–IV. Generally, deterministic equations for the echo magnitude [Eqs. (4) and (6)], which are based on solutions to the wave equation, are randomized with respect to various physical quantities such as random location in beampattern and random orientation of scatterer. They are randomized using fundamental statistics equations given in Sec. IV C. This approach differs from other approaches such as first randomizing parameters of the governing differential equation before solving the equation [see, for example, the summary in Sec. 12.6 of Jakeman and Ridley (2006), and references therein].

In Sec. VI, the treatment begins with the simplest of cases—single-frequency signals of infinite extent in which echoes from all scatterers completely overlap and direct path propagation in a homogeneous medium without the complicating effects due to reflections from boundaries. These idealized signals approximate gated sine waves where the echoes overlap significantly. In addition, given the importance of the random modulation effect of the beampattern on the echo when a scatterer is randomly located in the beam, the echo statistics are first described without those effects. This is equivalent to the sensor system having an omnidirectional beam. In Sec. VII, beampattern effects are then incorporated into the analysis. In Sec. VIII, the work is further extended into more realistic and complex cases involving pulsed signals (narrowband and broadband) in which echoes may only partially overlap, the presence of a boundary near a scatterer, and propagation and scattering in a waveguide with a heterogeneous medium.

The number of combinations of types of systems, signals, signal processing and beamforming algorithms, and environments is limitless and cannot be adequately described within this tutorial. The material is therefore aimed toward the more fundamental aspects of echo statistics that are not specific to any particular system or environment, but that many applications either have in common or could use as a basis. For example, most material involves the following: (1) Signal type: Long, narrowband signals are used in which the echoes from all scatterers completely overlap (in all cases in Secs. VI and VII, and some cases in Sec. VIII). (2) Signal processing: The magnitude of the signal as measured at a single instant in time is measured—that is, “first-order statistics” is modeled (in all cases in Secs. VI–VIII). The instant in time may be fixed or randomly selected, but it is not adaptively chosen according to a particular echo magnitude. (3) Processing of beam data and/or beamforming: Echoes from a single beam are modeled (fixed or scanning; or one selected from a multi-beam system) in which the scatterers are randomly distributed in space (in all cases in Secs. VII and VIII). The echo is sampled from a single beam for a random spatial distribution of scatterers and the beam is not steered adaptively to select or focus on a particular scatterer. (4) Environment: direct path geometries in which the medium is homogeneous and there are no reflecting boundaries (in Secs. VI and VII; and one case in Sec. VIII).

The following cases for systems, signals, and environments of greater complexity are examined in Sec. VIII: (1) Signal type and signal processing: pulsed signals are modeled in one example in which the echoes from the scatterers only partially overlap (Sec. VIII A). In Sec. VIII A, the pulse is further shortened through use of matched filter processing. (2) Environment: in two examples, geometries in which there are one or more boundaries present and there are heterogeneities in the medium (Secs. VIII B and VIII C, which cover boundary interaction and waveguide effects).

The several cases modeled in Sec. VIII, while far from spanning the many possible complex scenarios, provide examples for how the fundamental formulations involving the more simple cases can be applied to the more complex cases. For example, the case involving a pulsed signal shows how a time series can be constructed due to the interference between the partially overlapping echoes from the scatterers (Sec. VIII A). The examples involving the presence of one or more boundaries and/or heterogeneities in the medium show the different types of effects associated with the boundaries and medium heterogeneities. For a single boundary near a scatterer, it can provide an added source of interference due to the interaction of the incident signal and the boundary and scatterer (Sec. VIII B). For two parallel boundaries and/or a medium with a local minimum in wave speed, a waveguide is formed and the signal can propagate along multiple paths that are guided by the boundaries or local minimum (Sec. VIII C). These multiple paths represent additional sources of fluctuations in the echoes.

There are many more important cases not covered explicitly in this tutorial, but that can be described using these formulations. For example, there are systems that use peak sampling signal processing, such as recording the maximum echo magnitude in a time window or adaptively steering a beam toward the scatterer with the largest echo. When multiple scatterers are present, the process of peak sampling in both cases (time and angle/space) will bias the statistics toward higher values than the magnitudes modeled with first-order statistics and a fixed beam. This process involves “extremal” statistics (Stanton, 1985), which is outside the scope of this tutorial. However, for the time-based peak sampling, the method in Sec. VIII A that produces a time series could then serve as a basis of the time series with which a peak sampling algorithm could be applied. As shown in Stanton (1985), the bias in this case increases with the ratio of window duration to ping duration of the signal. For the case of a scanning beam or multi-beam system adaptively focusing on the peak echo in a field of multiple scatterers, this extremal statistics formulation can be adapted from a time series to a space series of echoes scanned across angles. The bias here will increase as the angular window is increased.

Another important case involves use of pulsed signals in the presence of boundaries in which the echoes from the scatterers are only partially overlapping. Here, the method given in Sec. VIII A to produce a time series can be incorporated into the formulations in Secs. VIII B or VIII C that involve the presence of boundaries. Furthermore, advanced signal processing such as peak sampling can also be incorporated into this case as described above.

Given the complex effects of the beampattern on the echo statistics (Fig. 4), fluctuations of the echoes without the influence of the beampattern are first examined separately. This is equivalent to a sensor system with an omnidirectional beam so that the echo value is the same regardless of angular location (at constant range) in the beam. In this simplified case, fluctuations can be due to a combination of interference between overlapping echoes and scattering effects. Fluctuations due to those effects are treated separately below.

This initial treatment solely involves analyses of generic random signals, devoid of explicit representation of the sensor system and field of scatterers. This analysis forms the foundation for the more complex analyses given later that first involve scattering from objects and, eventually, system effects. The case of an arbitrary number of arbitrary signals is first presented, which is then followed by some commonly used special cases, including the Rice and K PDFs.

When N arbitrary complex signals completely overlap, the resultant signal A is the coherent sum

where $ai$ and $Δi$ are the amplitude and phase of each individual signal, respectively, and are both considered arbitrary random variables. Since this equation models sinusoidal signals, all with the same frequency, the term $e−jωt$ that each signal has in common, has been suppressed as in the previous formulations.

Since $ai$ and $Δi$ are random variables, then so is A. The fluctuations of A from realization to realization depend strongly on the statistical properties of $ai$ and $Δi$⁠. The phase shifts $Δi$ play a major role in the fluctuations. For example, for the simple case in which $ai$ is constant for all i (i.e., $ai=a$⁠) and $Δi$ is randomly and uniformly distributed over the range 0–2$π$⁠, A will fluctuate greatly from realization to realization due to variability in constructive and destructive interference effects associated with phase variability alone (Figs. 2 and 6). In one realization, there may be complete constructive interference and A is at a maximum. In another realization, there may be complete destructive interference and A is at a minimum. And, generally, A will take on intermediate values due to partial interference.

The characteristics of the fluctuations in this case also depend greatly on N as illustrated. When there is only one signal (N = 1), the signal is single valued and the PDF of the signal magnitude $|A|$ is the delta function (Fig. 2). In the other extreme in which there are an infinite number of signals (N = $∞$⁠), the PDF of $|A|$ is the Rayleigh PDF as given in Eq. (21) (Fig. 2). The PDF of $|A|$ takes on other shapes for intermediate values of N (Fig. 6). In the more general case in which $ai$ is a random variable (not equal), the curves will fluctuate in a similar fashion and, in the limit of N = $∞$⁠, the PDF of $|A|$ will become Rayleigh via the CLT.

Equation (24) is broadly applicable to the scattering problem as it could represent the summation of scattering highlights from within a single scatterer or the summation of echoes from multiple scatterers. Characteristics of the scatterer, sensor system, and scattering geometry can be incorporated into N, $ai$⁠, and $Δi$⁠.

Equation (24) can be manipulated to model the important case of a signal in the presence of noise. In this case, one of the amplitudes $ai$ is held fixed while the others are randomly varied. Pulling the fixed-amplitude signal out of the summation, Eq. (24) is rewritten,

where the second term on the right-hand-side, representing noise, is the sum of a large number of random amplitude, random phase signals. The PDFs of the magnitudes of the fixed amplitude and summed signals are the delta function and Rayleigh PDF, respectively. Here, N must be sufficiently large so that the magnitude of the summation converges to a Rayleigh random variable. The PDF of $|A|$ is the Rice PDF (Rice, 1954),

where

and

The term $γ$ is the ratio of the mean squared values of the sine wave and noise (i.e., the power signal-to-noise ratio or SNR where the “signal” is the sine wave in this context), $σn$ is the rms value of the noise term, $⟨|A|2⟩$ is the mean square of the sine wave plus noise [i.e., mean square of Eq. (25)], and $I0$ is the zeroth-order modified Bessel function of the first kind. Note that in the original derivation by Rice, the noise term in Eq. (25) can also involve the summation of signals of different frequencies as well (i.e., more general than the single frequency case shown here).

The shape of the Rice PDF depends strongly on $γ$ (Fig. 7). For example, in the limit as $γ$ approaches infinity, the PDF is close to a Gaussian PDF. This corresponds to the limit of high SNR where the signal is dominated by the constant sine wave. In the other extreme, as $γ$ approaches zero, the PDF approaches the Rayleigh PDF. This latter case corresponds to the limit of low SNR where the signal is dominated by the noise. The shape of the Rice PDF changes smoothly for all intermediate values of $γ$⁠.

In addition to this formula being widely applicable in modeling noisy signals, it can also be applied to scattering problems where the “signal” [sine wave in Eq. (25)] is the mean scattered field and the “noise” [summation term in Eq. (25)] is the component of the scattered signal that fluctuates about the mean (Stanton and Clay, 1986; Stanton and Chu, 1992). For example, it could be used to model the fluctuations of the echo from a sphere near a rough interface where the individual echoes from the sphere and rough interface are delta-function- and Rayleigh-distributed, respectively. Or, rather than a rough interface, the sphere could be surrounded by a cloud of smaller scatterers whose individual echoes are of random phase. If the constant sine wave in Eq. (25) is considered to represent a single scatterer of interest (such as the sphere) and the noise term represents the background reverberation (such as from the rough interface or cloud of smaller scatterers), then the term $γ$ is the signal-to-reverberation ratio.

A more complex, but commonly occurring, case is when N and/or $ai$ in Eq. (24) are random variables. This can be divided into two categories—one in which $Δi$ is randomly and uniformly distributed over the range 0–2$π$⁠, and the other in which $Δi$ is non-uniformly distributed. The former category will first be discussed. If N follows a negative binomial PDF and its average value tends to infinity, the statistics of $|A|$ in Eq. (24) for arbitrary $ai$ is described by the K PDF (Jakeman and Pusey, 1976; Abraham and Lyons, 2002b),

where K is the modified Bessel function of the second kind (and served in the naming of the PDF), Γ is the gamma function, $αK$ is the shape parameter, and $λK$ is a scale parameter equal to the mean square of the signal divided by $αK$⁠.

The K PDF has a single mode and varies smoothly with $|A|$ (Fig. 8). The distribution tends to the Rayleigh PDF in the limit of high $αK$⁠. Note that the function K also has other, less common names, including the Basset function and the modified Bessel function of the third kind.

The K PDF can also be derived in several other ways. For example, it has been shown that for finite N and if $ai$ follows an exponential PDF, then $|A|$ is K-distributed (Abraham and Lyons, 2002b). Beyond methods involving summing sinusoidal signals, the K PDF has been shown to be due to the product of two independent random variables: a Rayleigh-distributed term and one that is chi distributed (Ward, 1981). This product was later written in extended form as the product of a Rayleigh-distributed term and the square root of a term that is gamma distributed (Abraham and Lyons, 2002b). Here, the square root of the gamma-distributed term is related to the chi-distributed term through analytical continuation of the integer number of summed terms in the chi distribution to a non-integer number in the gamma distribution. In another derivation, the K PDF has also been shown to result from a Rayleigh PDF whose mean-square value is gamma distributed (Jakeman and Tough, 1987). Both of these latter derivations are referred to as a “compound representation.”

Equation (29) has been widely used to describe echo statistics in both acoustic and electromagnetic applications. While there has generally not been a direct connection to the physics of the scattering, the tail of the distribution has generally followed those from experimental data after an empirical fit. As discussed above, through the interpretation of Abraham and Lyons (2002b), the number of scatterers has been related to the shape parameter of the K PDF for the specific case in which the amplitudes $ai$ of the individual echoes are exponentially distributed. Since $ai$ in this case is observed through the receiver of the sensor system, the exponential PDF includes the effects of both fluctuations from the stochastic nature of the scatterer and the variability due to the scatterer being randomly located in the beam. Details of those effects are given in Sec. VII. Also, in his expression of a K-distributed magnitude being due to the product of Rayleigh- and chi-distributed random variables, Ward (1981) attributed the Rayleigh term as being due to quickly varying interference between scatterers (such as from phase shift differences within a patch of scatterers) and the chi term being due to slowly varying changes in the echo from larger-scale variations in the “bunching” or patchiness of scatterers.

The above K PDF involves signals who phases are uniformly and randomly distributed over [0 2$π$]. However, when the distribution of phases is non-uniform, then $|A|$ can be described by the generalized K-distribution (not shown) (Jakeman and Tough, 1987). This distribution is described by three parameters—$λK$ and $αK$ as given above, plus a third that describes the non-uniform phase distribution. In addition, the generalized K PDF in Jakeman and Tough (1987) can also describe an n-dimensional random walk (generalized from the two-dimensional walk for the standard K PDF). For the case of a non-uniform phase distribution, this is a random walk with directional bias. Similar to the compound representation of the K PDF above being a Rayleigh PDF with its mean-square value being gamma distributed, with this directional bias, the generalized K PDF can be (compound) represented by a Rice PDF with the mean square noise and constant amplitude signal components each being gamma distributed in a correlated way (Jakeman and Tough, 1987). The generalized K PDF can be applicable to the case in which one or several scatterers dominate the scattering from within a field of many scatterers, thus skewing the distribution of phases into one that is non-uniform (Ferrara et al., 2011).

Another generalization of the K PDF is the homodyned K-distribution (not shown). Like the generalized K PDF, it is also a three-parameter distribution. And, like the generalized K PDF, it can also be (compound) represented by a Rice PDF, but with only the mean square noise component (but not the constant amplitude component) being gamma distributed (Jakeman and Tough, 1987).

The properties and potential relations to scattering of the K PDF, generalizations of the K PDF, and other generic PDFs are summarized in Destrempes and Cloutier (2010).

As discussed in Secs. IV C 6 and VI A 1, in the limiting case of the sum of an infinite number of random complex variables, the PDF of the magnitude of the sum is Rayleigh distributed (as per the CLT). Specifically, for Eq. (24) where the phase shifts $Δi$ are randomly and uniformly distributed over the range 0–2$π$⁠, in the limit of N = $∞$⁠, $|A|$ is Rayleigh distributed. Now, consider the case in which there is an ensemble of statistically independent realizations of the complex signal A (in the limit of large N) where the magnitude of each realization of A is Rayleigh distributed and with the same mean square value. Then by extension of the CLT, the magnitude of the sum of those realizations of A is likewise Rayleigh distributed, even for a finite number of realizations. However, if each realization of A is modulated by a multiplicative term with a magnitude that is a random variable, the resultant magnitude of the sum of a finite number of these (complex) signals can possibly be strongly non-Rayleigh. As described in Sec. VI A 3, such a modulation can be caused by patchiness of the scatterers, resulting in echo PDFs that can be derived through a compound representation. Also, when the scatterers are uniformly distributed (i.e., no patchiness effects), there can also be modulation caused by the beampattern where the echo from the scatterer is randomly modulated by its random location in the beam, causing strongly non-Rayleigh echoes as described in Secs. VII B and VII C.

Scatterers can range in complexity from the simplest of form, a point scatterer, in which the scattering amplitude is constant for all orientations, to an arbitrarily shaped object whose echo varies from orientation to orientation. In this section, the statistics of scattering by an individual are examined in a progression of complexity. The point scatterer, Rayleigh scatterer (defined below), and smooth and rough prolate spheroids with both fixed and random orientation are modeled. A summary of scatterers with other complexities are given at the end.

The simplest case is the point scatterer or, more precisely, very small scatterer. The dimensions of this scatterer are sufficiently small compared with an acoustic or electromagnetic wavelength that the echo is constant for all orientations. It is constant because there is only one scattering highlight from this object and nothing else with which to interfere. The PDF of the echo magnitude is the delta function [Figs. 2(a) and 3(a)].

At the other extreme from a point scatterer is a scatterer whose echo is Rayleigh distributed. This so-called “Rayleigh scatterer” can be in many forms. For example, it could be a small patch of many point scatterers, each with an echo whose phase is randomly and uniformly distributed over the range [0 2$π$]. Or, it could be a single spherical scatterer whose surface is randomly rough and can be described as a collection of many scattering highlights bounded by the surface. The echo from each highlight has a phase that is randomly and uniformly distributed over the range [0 2$π$]. In each case, from the central limit theorem, the echo from the patch or rough sphere is Rayleigh distributed [Fig. 2(c); Eq. (21)]. The fluctuations occur from ping to ping as the patch or rough sphere are rotated, or from realization to realization of a randomized spatial distribution of scatterers in the patch or randomized roughness of the rough sphere.

A sequence of formulations is presented, beginning with the deterministic description of the scattering by a smooth prolate spheroid at fixed orientation, then randomizing its orientation, and then further randomizing the spheroid by roughening its boundary. In contrast to the Rayleigh scatterer, the echoes from the randomized prolate spheroid are generally non-Rayleigh because of the elongated shape of the prolate spheroid. The degree to which the echoes are non-Rayleigh can be connected to various parameters of the scatterer through the physics-based formulas given below. The formulas given below are adapted from Bhatia et al. (2015).

We begin with a deterministic description of the scattering by a smooth prolate spheroid at fixed orientation. The spheroid is modeled as being impenetrable (acoustically “hard” or “soft” or, with an electromagnetic signal, perfectly conducting). Also, the scattering is in the “geometric optics” or high frequency limit where the acoustic or electromagnetic wavelength is much smaller than any dimension of the spheroid. For simplicity, only echoes from the front interface are analyzed and waves that travel around the boundary (i.e., circumferential waves) are ignored. Using the Kirchhoff approximation and the stationary phase approximation, the magnitude of the backscattering amplitude of the spheroid is (Chap. 4 of Crispin and Siegel, 1968)

where the “SS” subscript to the scattering amplitude refers to “smooth spheroid” and β is the angle between the direction of the incident acoustic or electromagnetic wave and the plane that is normal to the lengthwise axis of the prolate spheroid [Fig. 9(a)]. Note that this equation is equivalent to Eq. (7) of Bhatia et al. (2015), but is in a form in which no terms have singularities. β = 0 and β = $π/2$ correspond to broadside and end-on incidence, respectively, relative to the incident wave. The terms $c$ and $b1$ are the lengths of the semimajor and semiminor axes of the prolate spheroid, respectively (the length and width of the spheroid are $2c$ and $2b1$⁠, respectively). The spheroid is axisymmetric about the length-wise axis, leading to only one unique value of semiminor axis. The term $b1$ is not to be confused with the beampattern b. The aspect ratio of the spheroid is defined to be the ratio $c/b1$ (or, equivalently, length/width). The scattering is a strong function of the aspect ratio and orientation [Fig. 9(b)]. At broadside incidence, the above formula reduces to

We now randomize the scattering by first randomizing the orientation. This is done by making β a random variable with an associated PDF $pβ(β)$⁠. For simplicity, the prolate spheroid will only rotate in a single plane about its minor axis in a plane containing the direction vector of the incident plane wave [Fig. 9(a)].

Since the scattering amplitude $fss$ is a function of the random variable β, then $fss$ is a random variable as well. Inserting $|fSS|$ and $pβ(β)$ into Eq. (16), the PDF of the magnitude of the scattering amplitude of a randomly oriented smooth prolate spheroid is

Equation (16) was used because the scattering amplitude in Eq. (30) varies monotonically over the entire range of orientations. If the scattering amplitude varied non-monotonically, then Eq. (17) would have been required to calculate the echo PDF.

Inserting Eq. (30) into Eq. (32) gives the PDF of the magnitude of the scattering amplitude of a randomly oriented smooth prolate spheroid explicitly in terms of the dimensions of the prolate spheroid,

Calculation of the echo PDF requires knowledge of the orientation distribution. For the simple case in which the angles of rotation are uniformly and randomly distributed over the range [0 2$π$], the PDF of β is

where β only varies over the range 0–$π$/2 because of the symmetry of the scattering over the other angles.

The echo PDF for the randomly oriented smooth prolate spheroid is calculated for a range of aspect ratios (not including the case of a sphere where the aspect ratio is unity) by inserting Eq. (34) into Eq. (33) [Fig. 10(a)]. For each (non-unity) aspect ratio, the echo PDF is characterized by a smoothly varying function for most magnitudes, but with strong narrow peaks at the end points. These peaks are attributed to the fact that the backscattering near broadside and end-on incidence varies slowly with orientation angle, which increases the probability of occurrence at those corresponding echo values. The range of echo values as well as probability of occurrence are both shown to be a strong function of aspect ratio of the spheroid. For the case of a sphere, the echo is constant for all β and the echo PDF is a delta function [Fig. 10(a)].

We further randomize the scattering by roughening the boundary of the prolate spheroid, where the roughness here is the deviation from the mean boundary. In this case, the boundary is randomly rough throughout the entire surface of the spheroid. Furthermore, it is assumed to be sufficiently rough compared with an acoustic or electromagnetic wavelength such that for any orientation, the magnitude of the scattering amplitude is Rayleigh distributed. For the randomly rough surface, the echoes are assumed to be independent of each other from orientation to orientation (or realization to realization). In this limiting case, we model the scattering by being equal to the product of the magnitude of the scattering amplitude $|fss|$ of the smooth prolate spheroid at a particular orientation and an independent “modulation” random variable that follows a Rayleigh PDF. This latter modulation term has a unity rms and, for each orientation, the term is statistically independent from its value at all other orientations. For a randomly oriented prolate spheroid, $|fss|$ is also a random variable. Thus, the magnitude of the scattering amplitude for the randomly rough, randomly oriented prolate spheroid is the product of two random variables, $|fss|$ and the Rayleigh-distributed modulation term. The statistics of the echoes from randomly rough, randomly oriented prolate spheroids can be described using Eq. (20),

where the subscript “rs” refers to rough (randomly oriented) prolate spheroid and $pRay$ is the Rayleigh PDF of the modulation term [Eq. (21)]. The terms $|fss|$_min and $|fss|$_max are the minimum and maximum values of the magnitude of the scattering amplitude of the smooth prolate spheroid, respectively, which correspond to end-on and broadside incidence. Those two terms replaced the limits −$∞$ and +$∞$⁠, respectively, in the integral in Eq. (20) to reflect the fact that the magnitude of the scattering amplitude of the smooth prolate spheroid is within the range [$|fss|$_min, $|fss|$_max] and that the corresponding PDF $pss$ is zero for values of its argument outside that range. Without the above constraint, the integral still would have been constrained by the Rayleigh PDF whose argument is limited to only positive values, which would have led to integral limits of zero to infinity.

The resultant echo PDFs of the scattering by the randomly rough, randomly oriented prolate spheroid are significantly different from the smooth counterpart [Fig. 10(b)]. The curves for each aspect ratio are now smoothly varying, have a mode, and do not have any singularities.

There are many more important complexities of the scattering and associated echo statistics than described above and they are summarized briefly below. The above example involved the simplified case in which the boundary of the prolate spheroid was impenetrable, only the echo from the front interface was accounted for, and its roughness was sufficiently large so that the echo was Rayleigh distributed for any fixed orientation. However, it is not uncommon for an object to have circumferential and internally refracted/reflected waves and/or a boundary roughness that spans small-to-intermediate levels relative to a wavelength of the incident signal. Each of these latter cases can lead to statistics of the echo that are significantly different than that of the very rough impenetrable prolate spheroid discussed above and must be accounted for.

Given the complexity of predicting the scattering by objects with those various types of waves, the echo statistics will be summarized below according to a progression of complexity. In Sec. VI B 4 a, effects on the scattering due to circumferential and internal waves for smooth objects are briefly reviewed. In Sec. VI B 4 b, the object is roughened, and the statistics of the echoes from only the front interface are considered. This is similar to the case of the rough impenetrable prolate spheroid above but with the added complexity of accounting for small-to-intermediate roughness where the echoes are non-Rayleigh. In Sec. VI B 4 c, circumferential/internal waves associated with penetrable materials are incorporated with boundary roughness included.

Objects with common material properties not only produce a significant echo from the front interface, but they also support various classes of circumferential and internally refracted/reflected waves due to incident signals (acoustic: Hackman, 1993; Marston, 1992; Numrich and Uberall, 1992; Uberall, 1973; electromagnetic: Chen, 1998; Crispin and Siegel, 1968; Moser et al., 1981; Murphy et al., 1980; Newton, 1982; acoustic and electromagnetic: Bowman et al., 1987; Nussenzveig, 1969). For the case of impenetrable objects, such as infinitely dense materials for acoustics applications and perfectly conducting materials for electromagnetic applications, there is a weak signal that diffracts or “creeps” around the boundary and reradiates in all directions including back to the receiver. For the case of penetrable objects, such as elastic materials for acoustics applications and dielectric materials for electromagnetic applications, those diffracted waves also exist. In addition, with elastic objects, there can be strong surface elastic waves that are generated from the incident acoustic signal and will travel around the boundary (in addition to the diffraction boundary wave). Furthermore, other processes are also involved for the penetrable objects such as waves that refract into the interior of the object and reflect internally. Because of these various effects, the echo from a smooth symmetrical object such as a sphere or cylinder will generally be composed of the sum of echoes associated with the front-facing interface (sometimes called the “specular” echo), internal transmission and reflections, and circumferential waves.

The total echo from these objects will vary according to the constructive and destructive interference between the individual echoes associated with each of the different scattering phenomena described above. For a smooth object such as a sphere or cylinder, the phase of the echo associated with each phenomenon will vary strongly with signal frequency (or more precisely, ka, where a is the radius of the object). Because of these frequency dependences, which also vary with each phenomenon, the pattern of echo vs frequency will contain a series of peaks and strong nulls associated with the constructive and destructive interference, respectively, between the different phenomena.

The following discussion concerns the echo from the front interface only, in isolation from the circumferential and internal waves.

When the roughness of the surface is smaller than the wavelength of the incident signal or, more precisely, when k$σB$ < 1, where $σB$ is the standard deviation of the boundary (or, equivalently, the rms deviation from the mean boundary), the magnitude of the echo from the front-facing surface of a rough object at normal incidence to that surface is generally not Rayleigh distributed. This is because the phases of the echoes from the individual scattering features of the surface are relatively narrowly distributed, in contrast to being uniformly distributed [0 2$π$] such as in the above case of the (very) rough spheroid. This effect has been studied mostly in the context of scattering by rough planar interfaces and, to a much lesser extent, for rough bounded objects. In either case, in the limit of small roughness (i.e., k$σB$ ≪ 1), the echo PDF will tend to the delta function. In the opposite limit (i.e., k$σB$ ≫ 1), the echo PDF will tend to the Rayleigh PDF (or at least be Rayleigh-like) as in the above example.

It has been shown that the echo statistics from the full range of roughness at normal incidence of a planar boundary can be described by the Rice PDF, as given above in Eq. (26) (Stanton and Clay, 1986). In this formulation, the scattered signal is decomposed into the sum of two components, the coherent mean (constant amplitude) and fluctuation component. The PDFs of the magnitudes of the two components are the delta function and Rayleigh PDF, respectively. The mean component is related to the reflection coefficient of the interface modified by the term $e−2k2σB2$—a term originally derived for rough planar interfaces by Eckart (1953). The fluctuation term is related to the scattering phenomena that cause deviations in the echo from this mean. These two scattering terms, the constant amplitude and fluctuation components, correspond to the sinusoidal signal and noise, respectively, in the original Rice formulation.

The Rice PDF shape parameter has also been explicitly connected to parameters of the roughness and sensor system for rough planar interfaces (normal incidence), as summarized in Stanton and Clay (1986). Parameters of the roughness are the rms deviations (⁠$σB$⁠) from the normal to the mean surface and the two-dimensional autocorrelation function (along the surface) of the surface. Parameters of the sensor system are the frequency and beampattern. Although the studies summarized in Stanton and Clay (1986) are specific to underwater acoustic signals, the analysis involving the Rice PDF is formulated for scalar fields and could be applied to other types of sensor systems, such as radar, when scalar field representations are appropriate.

These formulations connecting the Rice PDF to fluctuations due to randomized scattering by a rough planar surface were extended to the case of acoustic scattering by rough elastic cylinders immersed in a fluid in Stanton and Chu (1992) and Gurley and Stanton (1993). Because of the complexity of the scattering by the elastic objects, the analysis in Stanton and Chu (1992) was divided into two formulations—a general one which described all echoes (i.e., including those associated with the front interface, the circumferential waves, and internally refracted waves) and a simple analysis involving only the front interface. The fluctuations of the echo from the rough front (curved) interface were related to the Rice PDF in the same manner as with the case of the rough planar interface involving a scalar field, but taking into account curvature of the surface. This simpler approach is clearly an approximation to the other, more general, case where the other waves played a major role in the fluctuations, especially near nulls due to interference effects (as discussed below). However, the Rice PDF, when used to model fluctuations of the echo from the front interface only, was demonstrated in these studies to predict the general trend of the fluctuations as parameters such as ka and k$σB$ were varied, where a is the mean radius of the rough cylinder.

An important element in the modeling of both the above simple and general cases was the fact that the Rice PDF assumes a noise term whose quadrature (i.e., real and imaginary) phase components have the same mean square values. However, the quadrature components of the random component of the scattered (scalar) signal from rough interfaces are generally not equal. Furthermore, the fluctuations of the scattered signal are sensitive only to the component of the “noise” (i.e., fluctuation component of scattered signal) that is in phase with the mean scattered field. These effects are quantified and accounted for in Stanton and Chu (1992) where the Rice PDF shape parameter is formulated in terms of the in-phase component of the fluctuation component of the scattered field.

The complexity is now increased by accounting for circumferential and internally refracting waves so that the echo consists of all waves—due to the front interface and circumferential/internal waves. As discussed above, in the case of the smooth penetrable object, these waves will interfere with each other and cause deep nulls in the pattern of echo magnitude versus frequency. When the penetrable object is roughened, the phases of the echoes associated with the various scattering phenomena will correspondingly vary, each in a different manner specific to the respective phenomenon. The total echo (sum of all components) will fluctuate from realization to realization due to the random roughness in a manner broadly similar to that from an impenetrable rough object but with important differences. For example, the center frequency of the nulls from the destructive interference will vary from realization to realization. Because of the steepness of the null (typically 30 dB variation within a narrow band of frequencies), the echo at frequencies near that of the null will fluctuate significantly.

The fluctuations near the null for acoustic scattering by randomly rough elastic cylinders immersed in a fluid have been observed in both numerical (Stanton and Chu, 1992) and experimental (Gurley and Stanton, 1993) studies. In these studies, the shape parameter of the Rice PDF, when fit with the simulations or experimental data, is shown to vary dramatically at frequencies near each null. The shape parameter decreases by as much as two orders of magnitude near the null which corresponds to a similar increase in the degree to which the echo fluctuates.

The roughness not only affects the phase shifts of circumferential and internal waves, but also their magnitude. In one study, the dominant acoustic Lamb-wave (or “plate wave”) of a spherical elastic shell along all meridional paths was randomized due to the roughness. The variability in path length and, hence, phase of this circumferential wave was related to the roughness. The total Lamb wave echo summed from all paths was shown to be attenuated exponentially due to the decrease in coherence of the signal (Stanton et al., 1998).

Although the above two examples involved scattering of acoustic waves by elastic objects immersed in a fluid, the same principles apply to other sensor systems such as for medical ultrasound or radar applications. In general, when the phases of these different types of scattered signals vary randomly due to roughness, then the shape parameter of the echo PDF will vary significantly near any null in the echo versus frequency pattern. Similarly, any type of circumnavigating or internally refracting wave can experience attenuation due to the different ray paths adding incoherently.

Beampattern effects are now added to the above treatment which involved the statistics of echoes from scatterers in the absence of beampattern effects (i.e., equivalent to an omnidirectional beam). Once the scatterers are placed in a directional beam of the sensor system, the echo received by the system becomes modulated by the beampattern. If the location of the scatterer is random, then the modulation is correspondingly random. In this case, the beampattern function is a random variable with a PDF referred to as the beampattern PDF.

The effect of the beampattern on echo statistics can be profound, as discussed in Sec. II C and illustrated (Fig. 4), and accounting for it can be complex. For example, in the simplest case of a point scatterer whose scattering amplitude is delta function distributed, the resultant echo PDF due to the scatterer being randomly located in a directional beam has a trend that is approximately power-law with strong structure superimposed [Fig. 4(d)].

The below treatment begins with a general formula for the echo PDF due to a single scatterer randomly distributed in a beam. The properties of the scatterer, spatial distribution, and beampattern are all arbitrary in this formula. Following the general formula is a progression of examples beginning with the simplest of cases—a point scatterer randomly and uniformly distributed within only the mainlobe of the beam of an axisymmetric transducer (i.e., excluding sidelobes). The analysis then extends to complex scatterers and the entire beam (including sidelobes) and, finally, to the case involving a beampattern from an arbitrary transducer in which the beampattern is not axisymmetric and a spatial distribution of scatterers that is not uniform. After the treatment of single scatterers, the echo statistics associated with multiple scatterers in the beam is presented. In the first set of examples, all scatterers are identical, which is then followed by the more general case of assemblages of scatterers of varying types.

The formulations can be applied to single beam (fixed or scanning) and multi-beam systems, provided that the scatterer(s) are randomly located in the beam. The formulations are general and are not specific to any particular system. Although specific types of signal processing or adaptive beamforming that some systems incorporate are not modeled, the formulations presented herein can serve as a basis for the modeling as discussed in Secs. V A and V B.

For the case in which a single scatterer is randomly located in the beam of the sensor system at approximately constant range, its angular coordinates $(θ, ϕ)$ are random variables. Since the beampattern is a function of these random variables, the beampattern function $b(θ, ϕ)$ is also now a random variable (Sec. IV C 3) and can be described by the beampattern PDF, $pb(b)$⁠. Consider now a scatterer whose scattering properties are random (such as a randomly rough and/or randomly oriented elongated scatterer). The scattering amplitude is now a random variable and its magnitude can be described by the PDF $ps(|fbs|)$ (Secs. II B and VI B). The magnitude of the echo $ẽ$ received by the system is equal to the product of the magnitude of the scattering amplitude $|fbs|$ and the beampattern $b(θ, ϕ)$ [Eq. (4)]. With both of these latter two terms being random variables, then $ẽ$ is also a random variable (Sec. IV C 4). Using Eq. (20), the PDF of $ẽ$ can be written in terms of $pS$ and $pb$ as [Ehrenberg (1972)]

where now $|fbs|$ (=$ẽ$/b) is implicitly the argument of $pS$⁠.

Equations (36a) and (36b) are equivalent because of the commutative nature of the product of the two random variables in Eq. (20). The limits of ±$∞$ in the integral in Eq. (20) are reduced to the ranges [$ẽ$ $∞$] and [0 1] in Eqs. (36a) and (36b), respectively, since the values of the beampattern b only span the range [0 1] and the corresponding beampattern PDF $pb$ is zero outside that range. Finally, the range $0≤b≤1$ is used in each equation as they apply to the entire beampattern, whereas some applications later will involve only portions of the beampattern, such as for values of the mainlobe only above the highest sidelobe. In those cases, the limits in the integrals are modified accordingly.

These expressions, originally derived by Ehrenberg (1972) [using Eq. (20)] for echo intensity with identical form, are given for echo magnitude (i.e., not intensity) in Ehrenberg et al. (1981) and are also described in reviews in Stanton and Clay (1986) and Ehrenberg (1989). While use of one form over the other [Eq. (36a) vs Eq. (36b)] was not explained in the early papers, it is possible that one form may be more conducive for evaluation, such as in numerical integration (Bhatia et al., 2015).

The above integral relationship in Eqs. (36a) and (36b) between the echo PDF $pe(ẽ)$ and the PDFs of the magnitude of the scattering amplitude and beampattern function is completely general, as it applies to an arbitrary stochastic scatterer that is randomly located (at approximately constant range) in an arbitrary beampattern over an arbitrary spatial distribution. The constraint of the scatterers being at approximately constant range is consistent with Eq. (4) as these equations apply to the scatterers distributed within a thin shell at a nearly constant distance from the sensor system. This eliminates the range-dependent effects in analysis of the echo fluctuations.

The beampattern PDF depends not only on the beampattern function $b(θ, ϕ)$⁠, but also the PDF, $pθ,ϕ(θ ,ϕ),$ of the angular location of the scatterer in the beam. This probability, in combination with the beampattern function, determines the degree to which the echo is randomly modulated by the beampattern. For example, if the scatterer were fixed in the center of the beam, then $pθ,ϕ$ is a delta function peaked at $(θ ,ϕ)$ = (0, 0) and the echo is multiplied by unity for all realizations [Fig. 4(a)]. In the other extreme, if the scatterer were randomly located throughout the entire half-space, then $pθ,ϕ$ is finite for all $θ and ϕ.$ In this latter case, the echo is randomly modulated by all values of the beampattern resulting in a wide range of echo values, even for a point scatterer [Fig. 4(d)].

Two simple examples are treated here involving the scatterer being randomly and uniformly distributed in a half-plane (2D) and half-space (3D) at approximately a constant range. The 2D and 3D cases apply to geometries in which the sensor system is detecting scatterers that are distributed throughout a thin semicircular arc and a thin hemispherical shell, respectively. The 2D case may apply to geometries where (1) the transducer is a line with a beampattern that only varies in one plane, (2) the sensor system is located within a thin layer of scatterers and is looking along the layer from within that layer, or (3) the system is in a waveguide and long-range echoes only vary with respect to one dimension (such as in the plane parallel to the waveguide boundaries as described in Sec. VIII C). In all of these cases, only a half-plane or half-space are considered as the transducer is assumed to be baffled sufficiently so that there is no “back” radiation.

In each example, $pθ,ϕ$ (which reduces to $pθ$ in some cases) is calculated using Eq. (8), where $pX(x)$ in that equation is used to represent the probability of occurrence of the scatterer per unit volume at angular location $(θ ,ϕ)$ and x represents the volume V (Chap. 10 of Medwin and Clay, 1998). The term $dPX$ $(=dPV)$ in Eq. (8) is the differential probability of occurrence of the scatterer in the differential volume dV at location $(θ ,ϕ)$⁠. For the case of the scatterers being located within a thin hemispherical shell of constant radius, which is the 3D geometry in most examples in this tutorial, then dx = dV = r²$sin θdθdϕ$Δr, where r is the radius of the hemispherical shell, $θ$ is the spherical polar angle, $ϕ$ is the spherical azimuthal angle, and Δr is the thickness of the shell. For the case of the scatterers being located within a thin-shelled semi-circular arc of constant radius r, which is the 2D geometry sometimes used in this tutorial, then dx = dV = r$dθΔwΔr$⁠, where $Δw$ is the (thin) width of the circular arc (strip). The total volume of the thin shell in each of these two cases is 2$π$r²Δr and $π$r$Δw$Δr, respectively.

For a scatterer uniformly distributed within the volume, then the differential probability of occurrence $dPV$ per unit differential volume $dV$ is held constant (i.e., $dPV/dV$ = p_v = constant). Using that constraint, and the fact that the integral of $dPV$ over the total volume that the scatterer can occupy is unity [and, hence, $dPV/dV$ = p_v = (total volume)⁻¹], then p_v = (⁠$π$r$Δw$Δr)⁻¹ and p_v = (2$π$r²Δr)⁻¹ for the 2D and 3D cases, respectively. Through these changes in variables, Eq. (8) becomes $dPV$ = $pθdθ$ and $dPV$ = $pθ,ϕdθdϕ$ for these two cases, respectively, where expressions for $pθ$ and $pθ,ϕ$ are given below.

For the 2D case in which the scatterer is randomly and uniformly distributed in a half-plane at approximately constant range, the probability density function of the angular location of the scatterer in spherical coordinates is determined using the above approach,

There is no dependence upon $θ$ since, at approximately constant range, the scatterer is uniformly distributed within a thin arc of nearly constant radius in that plane. There is no dependence upon $ϕ$ as it is fixed in this geometry. Note that although the polar angle $θ$ is normally restricted to the range $0≤θ≤π/2$⁠, it is varied over the range $−π/2≤θ≤π/2$ for fixed $ϕ$⁠. For the case of an axisymmetric beam centered at $θ$ = 0 which is typically the major response axis (MRA) of the beam, the expression $pθ=2/π$ for $0≤θ≤π/2$ has been used to eliminate redundant calculations (Bhatia et al., 2015).

For the 3D case in which the scatterer is randomly and uniformly distributed in a half-space at approximately constant range, the probability density function of its angular location is determined using the above approach (Medwin and Clay, 1998),

Although the scatterers are located throughout all values of $ϕ$⁠, $pθ,ϕ$ still does not depend upon $ϕ$ (as in the 2D case above). However, now $pθ,ϕ$ depends upon $θ$ because in this 3D polar-spherical coordinate system, the scatterer is randomly and uniformly distributed within a thin hemispherical shell in the range 0 ≤ $θ$ ≤ $π$/2. Calculations in this coordinate system involve annular rings (at constant spherical radius), each located at some angle $θ$ with a width of d$θ$ and spanning all values of $ϕ.$ Since the scatterer is randomly and uniformly distributed across all values of $ϕ$ within each ring, then the probability $pθ,ϕ$ only depends upon the area of each ring, which is proportional to $sin θ$ (which appears in the expression for dV above). Accounting for the uniform distribution across all $ϕ$ [0, $2π$] for a given value of $θ$ yields a factor (⁠$2π$⁠)⁻¹ in Eq. (38).

For the case in which the beampattern is symmetrical about the $θ=0$ axis, $pθ,ϕ$ in Eq. (38) can be integrated over all $ϕ$ [0, $2π$] for the simplified result

Calculating the PDF of the beampattern function depends upon the complexity of the beam. The simplest case is first examined involving an axisymmetric beam (i.e., due to a circular planar transducer) in which only the portion of the mainlobe above the value of the highest sidelobe is used [Figs. 4(c) and 11]. This results in the beampattern only being dependent upon a single random variable (spherical polar angle, $θ)$ and furthermore varying monotonically with respect to this angle. Using only the values of the mainlobe above that of the highest sidelobe is not only easier to calculate, but it also generally relates to the highest values of the echoes, which have a higher chance of being detectable above the system noise levels. With these simplifications, Eq. (16) can be used to describe the beampattern PDF for the case in which the beampattern is monotonic and dependent upon only a single random variable,

where $bSL$ is the value of the highest sidelobe of the beampattern and the notation $pθ(θ)$ represents $pθ,ϕ$ for this case where the scatterer distribution does not depend upon $ϕ$⁠. This corresponds to the scenario where a scatterer is uniformly and randomly distributed either in 2D or 3D as shown in Sec. VII A 2, but with the restriction that $−θSL≤θ≤θSL$ (fixed $ϕ$⁠) and $0≤θ≤θSL$ (all $ϕ$⁠) for the 2D and 3D cases, respectively, where $θSL$ corresponds to $bSL$ [Figs. 4(c) and 11]. With this restriction, $pθ(θ)=1/(2θSL)$ and $pθ(θ)=sin θ/(1−cos θSL)$ for the 2D and 3D cases, respectively. These two latter expressions are calculated in a similar manner as for Eqs. (37) and (39), respectively, except that they involve use of a smaller volume, subtended by the angle $θSL$⁠, over which $dPV$ is integrated for normalization.

The beampattern function for a circular planar transducer is (Kinsler et al., 2000)

where $aT$ is the radius of the transducer and $J1$ is the Bessel function of the first kind of order 1. The square of the brackets corresponds to the fact that this is a composite, or two-way beampattern, being produced by the product of the transmit and receive beampatterns which are identical to each other.

Using $pθ(θ)=sin θ/(1−cos θSL)$ from above and Eq. (41) in Eq. (40), the beampattern PDF for an axisymmetric beam and associated with the values of the mainlobe above the highest sidelobe is

where $J2$ is the Bessel function of the first kind of order 2. Here, the scatterers are assumed to be uniformly distributed at approximately constant range within the mainlobe for polar angles in the range $0≤θ≤θSL$ (all $ϕ$⁠) (i.e., a 3D case).

The beampattern PDF, when plotted on a log-log scale, has been shown to have a negative and nearly constant slope for the higher values of beampattern (Ehrenberg, 1972) (Figs. 4 and 12 of this paper). This pattern corresponds to the portion of the mainlobe higher than the highest sidelobe and also occurs over a wide range of beamwidths (Ehrenberg, 1972). Under these conditions, the beampattern PDF can be approximated using the following equation (Ehrenberg et al., 1981):

where the normalization constant k₀ ensures the integral of Eq. (43) over b is unity and is given by

Using Eq. (42) in the limit of b approaching unity or $kaT sin θ$ approaching zero (that is, for angles near the center of the beam), the exponent in Eq. (43) is (Chu and Stanton, 2010)

which shows that the slope of the beampattern PDF on a log-log plot varies only with $kaT$ (related to beamwidth) and is independent of b under these limiting conditions. This equation was derived under the assumption that the scatterers are uniformly distributed at approximately constant range within the mainlobe for polar angles in the range $0≤θ≤θSL$ (all $ϕ$⁠) (i.e., a 3D case). Note that the power-law form in Eq. (43) applies to both an intensity-based analysis (Ehrenberg, 1972; Ehrenberg et al., 1981) and magnitude-based analysis such as in this tutorial and in Chu and Stanton (2010).

Once the entire beampattern is accounted for in the echo statistics, the sidelobes become a significant factor. In this case, there can be more than one angle at which the beampattern achieves a certain value [Figs. 4(d) and 11]. In this non-monotonic case, Eq. (17) is now used to calculate the beampattern PDF,

Now, the full range of polar angles is used [$−π/2≤θ≤π/2$ (fixed $ϕ$⁠) and $0≤θ≤π/2$ (all $ϕ$⁠) for the 2D and 3D cases, respectively]. The summand of this equation is the same as Eq. (40) which corresponds to the portion of the mainlobe above the highest sidelobe level and which is a monotonic function. With this summation over m in Eq. (46), each value of m corresponds to a portion of the beampattern that is monotonic and is associated with the structure of the mainlobe and sidelobes (Fig. 11).

For a scatterer being uniformly distributed in a hemispherical shell, Eq. (39) is used for its spatial distribution $pθ(θ)$⁠. For the case of a circular planar transducer, Eqs. (39) and (41) are inserted into Eq. (46) to obtain the beampattern PDF for the entire beampattern,

The summand of this equation is similar to Eq. (42), which describes the PDF of the portion of the mainlobe above the highest sidelobe, and differs only by the term in Eq. (42) containing $θSL$⁠. By setting $θSL$ = $π/2$ in Eq. (42), thus allowing the polar angle to vary over the entire range $0≤θ≤π/2$ (all $ϕ$⁠), Eq. (42) becomes identical to the summand in Eq. (47).

With the sidelobes accounted for in Eq. (47), the beampattern PDF has significant structure involving singularities (Fig. 12). Each sharp peak in the PDF is associated with the peak of a sidelobe, while the smoothly varying portion with a nearly constant slope at the higher values of b is associated with the portion of the mainlobe above the highest sidelobe as discussed above. The beampattern PDF is also shown to vary with beamwidth (i.e., different $kaT$⁠). The narrower the beam, the more sidelobes are present, which leads to correspondingly more structure in the PDF. There is also some similarity in the occurrence of the sharp peaks as beamwidth is varied.

For any distribution of scatterers containing at least most of the main lobe of the beampattern, the beampattern PDF will generally be qualitatively similar for all distributions of scatterers. Specifically, the PDF will generally have a downward trend, such as with the (approximately) power law illustrated in Figs. 4 and 12. Naturally, there will be some differences associated with the different distributions. As shown in Fig. 4, if the scatterers are only in the main lobe of the beam and do not encounter sidelobes, then the beampattern PDF decreases monotonically with amplitude. Once the sidelobes are encountered, then there is structure superimposed on this decreasing trend. Furthermore, the beampattern PDF will also vary depending upon whether or not the scatterers reside in a half-plane containing the MRA or the entire half space (Fig. 13). While the beampattern PDF is shown to be generally similar between the two cases, there is a singularity associated with the 2D case at the maximum value of amplitude.