[Robotics] Sonar Noise Analysis

Sonar noise have two independent characteristics:

1. Intensity noise: usually multiplicative Gamma speckle.
2. Range noise: usually additive, spatially correlated residuals.
3. Time varying gain : the receiver increases gain for later-arriving echoes because later arrival corresponds to greater range.

Below, section 1 - 7 details the rationale behind applying gamma speckle multiplicatively to individual beams. Section 8 is how to introduce spatial correlation between beams using a Gaussian Kernel. Section 9 is about how to model TVG (time-varying gain) at the receiver end

1 - Phasor

Euler’s formula of a single complex number

A sinusoidal echo can therefore be represented by the complex number

where

For example, with (X=3) and (Y=4),

The physical waveform is the real part of the rotating complex signal:, whose frequency is fixed \(\omega\):

Expanding the above:

So any sinusoid can be represented as the real part of a complex phasor, a compact representation.

2- Rayleigh Distribution

A Rayleigh distribution describes the magnitude of a 2D vector whose two components are independent zero-mean Gaussian random variables. A sonar echo can be represented as a phasor:

Assume X and Y have the same standard deviation (important assumption):

The measured echo envelope is the magnitude is:

That magnitude (A) follows a Rayleigh distribution. Suppose one sonar resolution cell gives these in-phase and quadrature values:

Notice that Gaussian components can be positive or negative, but the magnitude cannot be negative:

A follows the Rayleigh probability density:

where \(\sigma\) is the standard deviation of X and Y. with (X) and (Y) independent. Then the joint distribution of (X,Y) is (w.r.t each realization (x,y)):

Using

we get

For (\(\sigma=1\)):

Mean and standard deviation

For a Rayleigh random variable,

and:

For (\(\sigma=1\)):

Therefore, coefficient of variation (CV) doesn’t depend on \(\sigma\):

3 - Sonar Pings Are Phasors

Inside one sonar resolution cell, the return may contain echoes from many small scatterers: sand grains, small rocks, etc. Each scatterer returns a complex acoustic wave with a different phase. Imagine many tiny scatterers inside one sonar beam:

1
2
3
4
sand grain 1 → small echo with phase 20°
sand grain 2 → small echo with phase 170°
sand grain 3 → small echo with phase 290°
...

Their complex echoes add:

When there are many scatterers with roughly random phases, the real and imaginary parts of (z) become approximately Gaussian. The real and imaginary parts of the sum become approximately Gaussian:

For example, if a Rayleigh-distributed amplitude has \(CV \approx 0.523\),

It means: If the logged sonar value is the raw linear echo amplitude from one independent look at a rough surface, its standard deviation should be about (52.3%) of its mean.

4 - Power Density Follows The Gamma Distribution

Define intensity or power as

For a transformation of random variables,

Where the derivative is:

Substitute the Rayleigh PDF:

This is an exponential distribution with mean

Therefore,

4 - 1 Gamma Distribution

On the other hand, A Gamma-distributed random variable (G) with shape (k) and scale (\(\theta\)) has PDF

Here:

• (k>0) is the shape parameter;
• (\(\theta\)>0) is the scale parameter;
• (\(\Gamma(k)\)) is the Gamma function.

For positive integers,

The mean and variance are

The coefficient of variation is

4 - 2 The exponential distribution is Gamma with k=1

Set

The Gamma PDF becomes

Since

we obtain

That is exactly the exponential PDF.

5 - Average of multiple pings produce a Gamma distribution

Suppose the sonar takes (L) independent measurements of the same underlying return:

Each single-look intensity is exponential:

Their sum is

The sum of (L) independent Gamma variables with the same scale is

Its PDF is

So the chain is

Usually the sonar averages the measurements:

If

then

Its variance becomes

Therefore, average more measurements reduces the relative fluctuation.

At pixel location ((r,c)), let the underlying clean intensity be

A single-look measured intensity is modeled as

where the second parameter is the scale. Because shape (=1), this is an exponential distribution with

Now we want to look at the normalized intensity distribution:

Then it’s

6 - Example with Hypothetical Numbers

As we increase the number of measurements of the same seabed surface:

1
2
3
4
5
L = 1   Exponential: strongly skewed, very noisy
L = 2   Still skewed, but smoother
L = 4   More concentrated near 1
L = 8   Considerably smoother
L → ∞   Concentrates around 1

Therefore, if we measured \(CV \in [0.41, 0.49]\) corresponds approximately to

7 - Apply Gamma Distributed Noise Multiplicatively On the Same Beam

Gaussian additive noise would say

That implies the same absolute noise strength for dark and bright returns and can produce negative intensities.

Gamma speckle says:

Where S is the Gamma noise. This better represents speckle because:

• intensity stays positive;
• bright returns fluctuate more in absolute magnitude;
• the relative variability is approximately constant;
• averaging multiple exponential looks naturally produces Gamma statistics.

The simplest model is independent Gamma speckle:

1
2
3
4
5
6
7
speckle = np.random.gamma(
    shape=L,
    scale=1.0 / L,
    size=clean_intensity.shape,
)

observed_intensity = clean_intensity * speckle

8 - Use Gaussian Kernel To Add Spatial Correlation To Noise

Suppose you start with white noise independently at every cell.:

Then filter it:

where (K) is a Gaussian kernel.

After convolution, adjacent cells share many of the same input random values:

1
2
output cell A uses nearby white-noise cells
output cell B uses almost the same nearby cells

Therefore, (g[r,c]) and (g[r,c+1]) become correlated.

Of course, an anisotropic kernel such as

1
sigma=(0.82, 0.75)

can create different correlation strengths in the two grid directions.

Then we add directly additive range noise:

One characteristic here is sonar signals are “heavy-tailed”. A regular gaussian noise is: 0.03, -0.04, 0.02, 0.06, -0.05, but a heavy tailed one is 0.03, -0.04, 0.02, 2.50, 0.06. “Tail” refers to the shape of the probability distribution, instead of its physical location.

In code, to model a gaussian noise with heavy tail

1
2
3
4
5
6
7
8
9
10
11
12
13
base = gaussian_blur(
    np.random.randn(rows, cols),
    sigma=(0.82, 0.75),
)

amplitude = np.exp(
    beta * gaussian_blur(
        np.random.randn(rows, cols),
        sigma=(4.0, 4.0),
    )
)

noise = base * amplitude

9 - TVG Modelling

TVG means time-varying gain. For synthetic data, TVG variation should usually be modeled as a smooth range-dependent multiplier, not independent random noise per pixel.

Substituting gives

where (p) is a residual exponent. For example:

Then:

• at (10) m, multiplier (=1);
• at (20) m, multiplier is

This means farther returns remain dimmer because TVG has not completely compensated the range loss.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
import numpy as np


def apply_simple_tvg_residual(
    clean_intensity: np.ndarray,
    clean_range: np.ndarray,
    reference_range_m: float = 10.0,
    residual_exponent: float = 1.2,
) -> np.ndarray:
    safe_range = np.clip(clean_range, 0.5, None)

    gain = (
        reference_range_m / safe_range
    ) ** residual_exponent

    return clean_intensity * gain

If the real data suggests gain drift or automatic gain adjustment, evolve parameters smoothly:

For example:

1
2
3
4
5
6
7
8
9
10
11
def update_tvg_exponent(
    previous: float,
    mean: float = 1.2,
    persistence: float = 0.98,
    innovation_std: float = 0.02,
) -> float:
    return (
        persistence * previous
        + (1.0 - persistence) * mean
        + np.random.normal(0.0, innovation_std)
    )

• Previous
  [Robotics] Time Series Analysis Metrics
• Next
  [Robotics] 3D-shape-detection-from-point-clouds