ARF of Beamforming

[OUCyf]

Greetings,
I'm Fu, and very nice package! When i use the beamforming method, I modify some of your code to get the array response function, but the results look like very weird.

According to the book 2018 seismic ambient noise, I set the cross spectral density C matrix to be the identity matrix in the following, I don't it is correct or not?

Best,
Fu

  # number of stations
  N = coordinates.shape[0]

  # compute steering_vectors
  frequency = np.mean(freq_band)
  n_vel, n_inc, n_azi, steering_vectors = compute_steering_vectors(
      coordinates,
      reference_receiver,
      dest_epsg,
      frequency,
      velocity,
      inclination,
      azimuth,
  )

  # compute covariance matrix
  C = np.ones((N, N), dtype=np.complex128)

  # compute beamforming power
  P: np.ndarray = np.einsum(
      "sn, nk, sk->s",
      steering_vectors.conj(),
      C,
      steering_vectors,
      optimize=True,
  )

  # reshape beamforming power
  P_reshape = np.reshape(P, (n_inc, n_azi, n_vel))
  P = np.real(P_reshape)

and the function compute_steering_vectors:

def compute_steering_vectors(
    coordinates,
    reference_receiver,
    dest_epsg,
    frequency: float,
    velocity: Union[float, tuple] = 6000,
    inclination: Union[float, tuple] = (-90, 90, 1),
    azimuth: tuple = (0, 360, 1),
) -> None:
    r"""Precompute the steering vectors

    Compute the steering vectors for the specified parameter range. For parameters that are specified as a tuple,
    the grid search is performed over the range: (min_value, max_value, increment)

    Parameters
    ----------
        frequency : :obj:`float`
            Discrete frequency at which beamforming is performed
        velocity : :obj:`float` or :obj:`tuple`
            Specifies the velocity as a float (if known) or grid over which search is performed
        inclination : :obj:`tuple`
            Specifies inclination grid over which search is performed
        azimuth : :obj:`tuple`
            Specifies azimuth grid over which search is performed

    """
    # number of stations
    N = coordinates.shape[0]

    # compute inclination, azimuth and velocity vectors
    if isinstance(velocity, tuple):
        velocity_gridded = np.arange(
            velocity[0],
            velocity[1] + velocity[2],
            velocity[2],
        )
    else:
        velocity_gridded = np.array([velocity])

    if isinstance(inclination, tuple):
        inclination_gridded = np.radians(
            np.arange(
                inclination[0],
                inclination[1] + inclination[2],
                inclination[2],
            )
        )
    else:
        inclination_gridded = np.radians(np.array([inclination]))

    azimuth_gridded = np.radians(
        np.arange(azimuth[0], azimuth[1] + azimuth[2], azimuth[2])
    )
    n_vel = len(velocity_gridded)
    n_inc = len(inclination_gridded)
    n_azi = len(azimuth_gridded)

    # create grid
    inclination_gridded, azimuth_gridded, velocity_gridded = np.meshgrid(
        inclination_gridded, azimuth_gridded, velocity_gridded, indexing="ij"
    )

    # compute relative coordinates, and convert lat/lon/elev to y/x/z
    utm_x, utm_y = latlon_2_utm(
        coordinates[:, 0], coordinates[:, 1], dest_epsg=dest_epsg
    )
    coords_utm = np.column_stack((utm_x, utm_y, coordinates[:, 2]))
    coords = coords_utm - np.tile(coords_utm[reference_receiver, :], (N, 1))
    coords = np.asmatrix(coords)

    # compute wave vector and wave number
    wave_vector_x = (np.sin(inclination_gridded) * np.cos(azimuth_gridded)).ravel()
    wave_vector_y = (np.sin(inclination_gridded) * np.sin(azimuth_gridded)).ravel()
    wave_vector_z = (np.cos(inclination_gridded)).ravel()
    wave_vector_x, wave_vector_y, wave_vector_z = (
        np.asmatrix(wave_vector_x).T,
        np.asmatrix(wave_vector_y).T,
        np.asmatrix(wave_vector_z).T,
    )
    wave_number = (-2 * np.pi * frequency / velocity_gridded).ravel()
    wave_number = np.asmatrix(wave_number).T

    # compute steering vectors, steering_vectors = exp(i * k * np.dot(wave_vector * coords))
    steering_vectors: np.ndarray = np.exp(
        1j
        * np.multiply(
            np.tile(wave_number, (1, N)),
            (
                wave_vector_x * coords[:, 0].T
                + wave_vector_y * coords[:, 1].T
                + wave_vector_z * coords[:, 2].T
            ),
        )
    )

    # ensure that steering vectors are unit vectors
    steering_vectors = steering_vectors / np.sqrt(N)

    return n_vel, n_inc, n_azi, steering_vectors

[OUCyf]

Or it could be possible to add the arf function in the package? Hope it will be easy for you!

[solldavid]

Dear Fu,

Apologies for my belated reply. At the moment, there is no direct way to plot array response functions in TwistPy. However, there is a simple workaround that allows you to readily compute ARF's with the current version of the code.

What you should do is generate some reference data for a monochromatic plane-wave with a specific reference slowness (typically, the reference slowness vector is chosen to be equal to 0, corresponding to a plane wave at vertical incidence) and then perform beamforming on this reference data. The output should then directly correspond to the ARF.

Please find an example for a cross-shaped array with 0 reference slowness below:

import numpy as np
import matplotlib.pyplot as plt
from obspy import Trace, UTCDateTime
from twistpy.array_processing import BeamformingArray
from twistpy.convenience import to_obspy_stream


N = 5  # Number of sensors per line
s_ref = np.array([0, 0, 0])  # Reference slowness vector
f = 0.02  # Frequency of the signal (Hz)
t = np.arange(0, 1 / f, 1 / f / 1001)  # Time axis stretching over 1 period
x = 1e3 * np.linspace(-300, 300, N)  # x-coordinates of the sensors
y = 1e3 * np.linspace(-300, 300, N)  # y-coordinates of the sensors
x = np.concatenate((x, np.zeros(x.shape)))
y = np.concatenate((np.zeros(y.shape), y))
z = np.zeros((2 * N,))  # z-coordinates of the sensors (Assuming a flat topography)

coordinates = np.array([x, y, z]).T

# Generate a monochromatic plane wave propagating in the direction given by the reference slowness vector
data = np.sin(
    2 * np.pi * np.dot(np.array([x, y, z]).T, s_ref)
    - 2 * np.pi * f * np.reshape(t, (t.size, 1))
)
start_time = UTCDateTime(2022, 2, 1, 10, 00, 00, 0)
data_st = to_obspy_stream(data.T, start_time, t[1] - t[0])


# Beamforming
inclination = (
    0,
    90,
    1,
)  # Search space for the inclination in degrees (min_value, max_value, increment)
azimuth = (
    0,
    360,
    1,
)  # Search space for the back-azimuth in degrees (min_value, max_value, increment)
velocity = 1000.0  # Intra-array velocity in m/s
# The array response function will be plotted as a function of
# the horizontaal slowness which is given by s = sin(inclination)* (1 / velocity)
s_max = 1 / velocity  # Maximum slowness (s/m)


# is unknown and part of the search
number_of_sources = 1  # Specify the number of interfering sources that will be estimated in the time window
# (only relevant for MUSIC)

array = BeamformingArray(name="Test Array", coordinates=coordinates)
array.add_data(data_st)
array.compute_steering_vectors(
    frequency=f,
    inclination=inclination,
    azimuth=azimuth,
    intra_array_velocity=velocity,
)

freq_band = (0.01, 0.03)  # Frequency band of interest
P_BARTLETT = array.beamforming(
    method="BARTLETT",
    event_time=start_time + t[501],
    frequency_band=freq_band,
    window=1,
)
P_MVDR = array.beamforming(
    method="MVDR", event_time=start_time + t[501], frequency_band=freq_band, window=1
)

P_MUSIC = array.beamforming(
    method="MUSIC",
    event_time=start_time + t[501],
    frequency_band=freq_band,
    window=1,
    number_of_sources=1,
)

# Plot the ARF for the different methods
azi_plot = np.arange(azimuth[0], azimuth[1] + azimuth[2], azimuth[2])
inc_plot = np.arange(inclination[0], inclination[1] + inclination[2], inclination[2])

azi_plot, inc_plot = np.meshgrid(
    np.radians(azi_plot), np.sin(np.radians(inc_plot)) * (s_max)
)

fig_bf_polar, ax_bf_polar = plt.subplots(
    1, 3, sharex=True, sharey=True, figsize=(15, 6), subplot_kw=dict(polar=True)
)
for ax_p in ax_bf_polar:
    ax_p.set_theta_direction(-1)
    ax_p.set_theta_offset(np.pi / 2.0)
ax_bf_polar[0].pcolormesh(azi_plot, inc_plot, P_MUSIC.squeeze())
ax_bf_polar[0].set_xlabel("Horizontal Slowness")
ax_bf_polar[0].set_title("MUSIC ARF")

ax_bf_polar[1].pcolormesh(azi_plot, inc_plot, P_MVDR.squeeze())
ax_bf_polar[1].set_xlabel("Horizontal Slowness")
ax_bf_polar[1].set_title("MVDR (Capon) ARF")

ax_bf_polar[2].pcolormesh(azi_plot, inc_plot, P_BARTLETT.squeeze())
ax_bf_polar[2].set_xlabel("Horizontal Slowness")
ax_bf_polar[2].set_title("BARTLETT (Conventional) ARF")

plt.figure()
plt.plot(x, y, "vk")
plt.title("Array geometry")
plt.xlabel("X (m)")
plt.ylabel("Y (m)")
plt.grid()
plt.show()

I didn't have time to test this extensively, so it would be great if you could confirm that this is what you are looking for. At first glance, the plots I get look as expected.

I will keep this issue open and potentially include a simpler way to plot ARF's in the future.

Best regards,

David

[OUCyf]

Thanks David! I appreciate that a lot!