Speckle Statistics

This notebook will provide an empirical demonstration of speckle - how it originates, how it visually and statistically looks like, and some of the most common approaches to filter it.

Speckle is defined as a kind of noise that affects all radar images. Given the multiple scattering contributions originating from the various elementary objects present within a resolution cell, the resulting backscatter signal can be described as a random constructive and destructive interference of wavelets. As a consequence, speckle is the reason why a granular pattern normally affects SAR images, making it more challenging to interpret and analyze them.

Credits: ESRI

import json
from functools import partial

import holoviews as hv
import intake
import matplotlib.pyplot as plt
import numpy as np
from holoviews.streams import RangeXY
from scipy.ndimage import uniform_filter

hv.extension("bokeh")

Let’s make an example of a cornfield (with a typical backscattering value of about -10 dB). According to the following equation:

\[ \sigma^0 = \frac{1}{\text{area}} \sum_{n \in \text{area}} \sigma_n \]

We should ideally have a uniform discrete sigma naught \(\sigma^0\) value, given that the cornfield pixel is the only individual contributor.

However, since we already learned from the previous notebooks that a pixel’s ground size can be in the order of tens of meters (i.e., 10 meters for Sentinel-1), we can imagine that different distributed targets in the scene contribute to the global backscattered information.

Let´s replicate this behavior with an ideal uniform area constituted by 100 pixels and then by adding 30% of speckle.

ideal_backscatter = -10  # in dB, a typical value for cornfields
width = 12
size = (width, width)
ideal_data = np.full(size, ideal_backscatter)
ideal_data_linear = 10 ** (
    ideal_data / 10
)  # Convert dB to linear scale for speckle addition

speckle_fraction = 0.3
num_speckled_pixels = int(
    size[0] * size[1] * speckle_fraction
)  # Rayleigh speckle noise
speckled_indices = np.random.choice(
    width * width, num_speckled_pixels, replace=False
)  # random indices for speckle

# Initialize speckled data as the same as the ideal data
speckled_data_linear = ideal_data_linear.copy()

speckle_noise = np.random.gumbel(scale=1.0, size=num_speckled_pixels)
speckled_data_linear.ravel()[
    speckled_indices
] *= speckle_noise  # Add speckle to the selected pixels

ideal_data_dB = 10 * np.log10(ideal_data_linear)
speckled_data_dB = 10 * np.log10(speckled_data_linear)
plt.figure(figsize=(16, 10))

# Ideal data
plt.subplot(2, 2, 1)
plt.imshow(ideal_data_dB, cmap="gray", vmin=-20, vmax=0)
plt.title("Ideal Backscatter (Cornfield)")
plt.colorbar(label="Backscatter (dB)")

# Speckled data
plt.subplot(2, 2, 2)
plt.imshow(speckled_data_dB, cmap="gray", vmin=-20, vmax=0)
plt.title(f"Speckled Backscatter ({int(speckle_fraction * 100)}% of Pixels)")
plt.colorbar(label="Backscatter (dB)")

bins = 25
hist_ideal, bins_ideal = np.histogram(ideal_data_dB.ravel(), bins=bins, range=(-20, 0))
hist_speckled, bins_speckled = np.histogram(
    speckled_data_dB.ravel(), bins=bins, range=(-20, 0)
)
max_freq = max(
    hist_ideal.max(), hist_speckled.max()
)  # maximum frequency for normalization

# Histogram for ideal data
plt.subplot(2, 2, 3)
plt.hist(ideal_data_dB.ravel(), bins=bins, range=(-20, 0), color="gray", alpha=0.7)
plt.ylim(0, max_freq)
plt.title("Histogram of Ideal Backscatter")
plt.xlabel("Backscatter (dB)")
plt.ylabel("Frequency")

# Histogram for speckled data
plt.subplot(2, 2, 4)
plt.hist(speckled_data_dB.ravel(), bins=bins, range=(-20, 0), color="gray", alpha=0.7)
plt.ylim(0, max_freq)
plt.title(f"Histogram of Speckled Backscatter ({int(speckle_fraction * 100)}%)")
plt.xlabel("Backscatter (dB)")
plt.ylabel("Frequency")

plt.tight_layout()

/tmp/ipykernel_3148/3459812397.py:26: RuntimeWarning: invalid value encountered in log10
  speckled_data_dB = 10 * np.log10(speckled_data_linear)

../../../_images/696848bf16601b4fff6a8d466849f74740ca45465ab0632933dbd04e563b0cda.png

Figure 1: Synthetic data that emulates speckles in microwave backscattering

We can imagine that the second plot represents a real SAR acquisition over a cornfield, while the first plot represents an ideal uniform SAR image over a cornfield land (no speckle). The introduction of a simulated 30% speckle noise could be related to the presence of distributed scatterers of any sort present in the scene, which would cause a pixel-to-pixel variation in terms of intensity.

All the random contributions (such as the wind) would result in a different speckle pattern each time a SAR scene is acquired over the same area. Many subpixel contributors build up a complex scattered pattern in any SAR image, making it erroneous to rely on a single pixel intensity for making reliable image analysis. In order to enhance the degree of usability of a SAR image, several techniques have been put in place to mitigate speckle. We will now show two of the most common approaches: the temporal and the spatial filter.

Spatial filtering

We first introduce a common spatial filter. The Lee filter is an adaptive speckle filter. The filter works using a kernel window with a configurable size, which refers to the dimensions of the neighborhood over which the filter operates. The kernel slides across the data, applying the smoothing operation at each pixel position of the image. It follows three assumptions:

SAR speckle is modeled as a multiplicative noise - the brighter the area the noisier the data.
The noise and the signal are statistically independent of each other.
The sample mean and sample variance of a pixel is equal to its local mean and local variance.

This approach comes with some limitations: it reduces the spatial resolution of the SAR image.

Let’s build up a function for applying a Lee filter with a kernel window size of 7 (do not forget to switch back to linear units before doing this computation and to dB after it):

def lee_filter(raster, size=7):
    """
    Parameters:
    raster: ndarray
        2D array representing the noisy image (e.g., radar image with speckle)
    size: int
        Size of the kernel window for the filter (must be odd, default is 7)

    Returns:
    filtered_image (ndarray): The filtered image with reduced speckle noise
    """

    raster = np.nan_to_num(raster)
    raster = 10 ** (raster / 10)

    # Mean and variance over local kernel window
    mean_window = uniform_filter(raster, size=size)
    mean_sq_window = uniform_filter(raster**2, size=size)
    variance_window = mean_sq_window - mean_window**2

    # Noise variance estimation (this could also be set manually)
    overall_variance = np.var(raster)

    # Compute the Lee filter
    weights = variance_window / (variance_window + overall_variance)

    return 10 * np.log10(mean_window + weights * (raster - mean_window))

plot_variability(fused_ds, filter_fun_spatial=lee_filter)

Figure 3: Lake Neusiedl \(\sigma^0_E\) with a Lee filter applied.

Temporal filtering

Temporal filtering would involve taking the average of all previous (past) observations for each pixel. This approach comes with some limitations: it takes out the content-rich information tied to the temporal variability of backscatter.

def temporal_filter(raster):
    """
    Parameters:
    raster: ndarray
        3D array representing the noisy image over time
        (e.g., radar image with speckle)

    Returns:
    filtered_image (ndarray): The filtered image with reduced speckle noise
    """

    return raster.mean("time")

plot_variability(fused_ds, filter_fun_temporal=temporal_filter)

Figure 4: Lake Neusiedl \(\sigma^0_E\) with a temporal filter applied.

Let´s observe the histograms of the two plots. Especially in the region around the lake, it is clear that the distribution is now less dispersed and more centered around a central value.

Speckle Statistics

Lake Neusiedl data

Spatial filtering

Temporal filtering