Matched-Field Processing - Synthetic 3D case

[1]:

from itertools import product
import dask.array as da
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
import pylab as plt
import torch
from tqdm import tqdm
from scipy.signal import ricker

# specific das_ice function
import das_ice.io as di_io
import das_ice.signal.filter as di_filter
import das_ice.processes as di_p
import das_ice.plot as di_plt
import das_ice.mfp as di_mfp

# classic librairy
import torch
import xarray as xr
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from tqdm import trange
%matplotlib inline

[2]:

from dask.distributed import Client,LocalCluster

cluster = LocalCluster(n_workers=8)
client = Client(cluster)
client.dashboard_link

[2]:

'http://127.0.0.1:8787/status'

Defined captor position

In this exemple, 4 sensors are defined.

[3]:

nb_sensor=4
random_sensor=[]
for i in range(nb_sensor):
    random_sensor.append([np.random.randint(-70, 71),np.random.randint(-70, 71),np.random.randint(-70, 71)])

sensors=np.array(random_sensor)

Generate signal

This section shows how to genrate signal that are similar to DAS data. Two events are compute.

[4]:

nb_sources=2
random_sources=[]
for i in range(nb_sources):
    random_sources.append([np.random.randint(-70, 71),np.random.randint(-70, 71),np.random.randint(-70, 71)])
    if i==0:
        dart=di_mfp.artificial_sources_freq(sensors,random_sources[-1],2500,window_length=.1,sampling_rate=25000)
    else:
        dart2=di_mfp.artificial_sources_freq(sensors,random_sources[-1],2500,window_length=.1,sampling_rate=25000)
        dart2['time']=dart2['time']+dart['time'][-1]+dart['time'][1]
        dart=xr.concat([dart,3*dart2],'time')

plt.figure()
for i in range(dart.shape[0]):
    dart[i,:].plot(label=i)

plt.legend()
plt.grid()
../_images/Tutorials_MFP_synthetic_case_3D_6_0.png

Performed MFP

This MPF implementation if taken from the dask implementation from this git repository.

3D MFP

Example of 3D MFP.

[5]:

MFP3D=di_mfp.MFP_3D_series(dart,0.1,sensors,xrange=[-70,70],yrange=[-71,71],zrange=[-72,72],dx=1,dy=1,dz=1)

/home/chauvet/miniforge3/envs/das_ice/lib/python3.11/site-packages/distributed/client.py:3362: UserWarning: Sending large graph of size 1.57 GiB.
This may cause some slowdown.
Consider loading the data with Dask directly
 or using futures or delayed objects to embed the data into the graph without repetition.
See also https://docs.dask.org/en/stable/best-practices.html#load-data-with-dask for more information.
  warnings.warn(

Show results

The MFP is finding the best match for:

[10]:

i_event=0
find_sources=MFP3D[...,i_event].where(MFP3D[...,i_event] == MFP3D[...,i_event].max(), drop=True).coords
find_sources

[10]:

Coordinates:
  * velocity   (velocity) int64 8B 2500
    true_time  float64 8B 0.0
  * x          (x) int64 8B 16
  * y          (y) int64 8B 69
  * z          (z) int64 8B -9

The true sources was in:

[11]:

random_sources[i_event]

[11]:

[16, 69, -9]

[12]:

matching_index = np.where(MFP3D.z.values == find_sources['z'].values)[0]
MFP3D[...,matching_index,0,i_event].plot()
plt.scatter(random_sources[i_event][1],random_sources[i_event][0],s=80,marker='s',color='k',label='random_sources')
plt.scatter(find_sources['y'],find_sources['x'],s=35,color='y',label='found with MFP')
plt.legend()

[12]:

<matplotlib.legend.Legend at 0x75a4e7249c90>
../_images/Tutorials_MFP_synthetic_case_3D_14_1.png

3D plot of the MFP Bartlett with a slider to plot the different plane of the Bartlett field.

[9]:

fig=di_mfp.plot_MFP3D(MFP3D[...,0,i_event],sensors,true_source=random_sources[i_event])
fig.update_layout(
    scene_camera=dict(
        eye=dict(x=1.2, y=1.2, z=1.2)
    )
)
fig.show("png")
../_images/Tutorials_MFP_synthetic_case_3D_16_0.png