.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "content/user-guide/tutorials/12-seismic/plot_fwd_1_tomography_2D.py"
.. LINE NUMBERS ARE GIVEN BELOW.
.. only:: html
.. note::
:class: sphx-glr-download-link-note
:ref:`Go to the end <sphx_glr_download_content_user-guide_tutorials_12-seismic_plot_fwd_1_tomography_2D.py>`
to download the full example code.
.. rst-class:: sphx-glr-example-title
.. _sphx_glr_content_user-guide_tutorials_12-seismic_plot_fwd_1_tomography_2D.py:
Forward Simulation for Straight Ray Tomography in 2D
====================================================
Here we module *simpeg.seismic.straight_ray_tomography* to predict arrival
time data for a synthetic velocity/slowness model. In this tutorial, we focus on the following:
- How to define the survey
- How to define the forward simulation
- How to predict arrival time data
.. GENERATED FROM PYTHON SOURCE LINES 16-19
Import Modules
--------------
.. GENERATED FROM PYTHON SOURCE LINES 19-33
.. code-block:: Python
import os
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from discretize import TensorMesh
from simpeg import maps
from simpeg.seismic import straight_ray_tomography as tomo
from simpeg.utils import model_builder
save_file = False
.. GENERATED FROM PYTHON SOURCE LINES 34-42
Defining the Survey
-------------------
Here, we define survey that will be used for the forward simulation. The survey
consists of a horizontal line of point receivers at Y = 100 m and a horizontal
line of point sources at Y = -100 m. The shot by each source is measured by
all receivers.
.. GENERATED FROM PYTHON SOURCE LINES 42-65
.. code-block:: Python
# Define the locations for the sources and receivers.
x = np.linspace(-100, 100, 11)
y_receivers = 100 * np.ones(len(x))
y_sources = -100 * np.ones(len(x))
receiver_locations = np.c_[x, y_receivers]
source_locations = np.c_[x, y_sources]
# Define the list of receivers used by each source
receiver_list = [tomo.Rx(receiver_locations)]
# Define an empty list to store sources objects. Define each source and
# provide its corresponding receivers list
source_list = []
for ii in range(0, len(y_sources)):
source_list.append(
tomo.Src(location=source_locations[ii, :], receiver_list=receiver_list)
)
# Define they tomography survey
survey = tomo.Survey(source_list)
.. GENERATED FROM PYTHON SOURCE LINES 66-72
Defining a Tensor Mesh
----------------------
Here, we create the tensor mesh that will be used to predict arrival time
data.
.. GENERATED FROM PYTHON SOURCE LINES 72-80
.. code-block:: Python
dh = 10.0 # cell width
N = 21 # number of cells in X and Y direction
hx = [(dh, N)]
hy = [(dh, N)]
mesh = TensorMesh([hx, hy], "CC")
.. GENERATED FROM PYTHON SOURCE LINES 81-89
Model and Mapping on Tensor Mesh
--------------------------------
Here, we create the velocity model that will be used to predict the data. Since
the physical parameter for straight ray tomography is slowness, we must define
a mapping which converts velocity values to slowness values. The model
consists of a lower velocity block within a higher velocity background.
.. GENERATED FROM PYTHON SOURCE LINES 89-122
.. code-block:: Python
# Define velocity of each unit in m/s
background_velocity = 3000.0
block_velocity = 1500.0
# Define the model. Models in SimPEG are vector arrays.
model = background_velocity * np.ones(mesh.nC)
ind_block = model_builder.get_indices_block(np.r_[-50, 20], np.r_[50, -20], mesh.gridCC)
model[ind_block] = block_velocity
# Define a mapping from the model (velocity) to the slowness. If your model
# consists of slowness values, you can use *maps.IdentityMap*.
model_mapping = maps.ReciprocalMap()
# Plot Velocity Model
fig = plt.figure(figsize=(6, 5.5))
ax1 = fig.add_axes([0.15, 0.15, 0.65, 0.75])
mesh.plot_image(model, ax=ax1, grid=True, pcolor_opts={"cmap": "viridis"})
ax1.set_xlabel("x (m)")
ax1.set_ylabel("y (m)")
ax1.plot(x, y_sources, "ro") # source locations
ax1.plot(x, y_receivers, "ko") # receiver locations
ax2 = fig.add_axes([0.82, 0.15, 0.05, 0.75])
norm = mpl.colors.Normalize(vmin=np.min(model), vmax=np.max(model))
cbar = mpl.colorbar.ColorbarBase(
ax2, norm=norm, orientation="vertical", cmap=mpl.cm.viridis
)
cbar.set_label("$Velocity (m/s)$", rotation=270, labelpad=15, size=12)
.. image-sg:: /content/user-guide/tutorials/12-seismic/images/sphx_glr_plot_fwd_1_tomography_2D_001.png
:alt: plot fwd 1 tomography 2D
:srcset: /content/user-guide/tutorials/12-seismic/images/sphx_glr_plot_fwd_1_tomography_2D_001.png
:class: sphx-glr-single-img
.. GENERATED FROM PYTHON SOURCE LINES 123-129
Simulation: Arrival Time
------------------------
Here we demonstrate how to predict arrival time data for the 2D straight
ray tomography problem using the 2D Integral formulation.
.. GENERATED FROM PYTHON SOURCE LINES 129-137
.. code-block:: Python
# Define the forward simulation. To do this we need the mesh, the survey and
# the mapping from the model to the slowness values on the mesh.
simulation = tomo.Simulation(mesh, survey=survey, slownessMap=model_mapping)
# Compute predicted data for some model
dpred = simulation.dpred(model)
.. GENERATED FROM PYTHON SOURCE LINES 138-141
Plotting
-----------------------------------------------
.. GENERATED FROM PYTHON SOURCE LINES 141-161
.. code-block:: Python
n_source = len(source_list)
n_receiver = len(x)
dpred_plotting = dpred.reshape(n_receiver, n_source)
fig = plt.figure(figsize=(8, 5))
ax = fig.add_subplot(111)
obs_string = []
for ii in range(0, n_source):
ax.plot(x, dpred_plotting[:, ii])
obs_string.append("source {}".format(ii + 1))
ax.set_xlim(np.min(x), np.max(x))
ax.set_xlabel("x (m)")
ax.set_ylabel("arrival time (s)")
ax.set_title("Positions vs. Arrival Time")
ax.legend(obs_string, loc="upper right")
.. image-sg:: /content/user-guide/tutorials/12-seismic/images/sphx_glr_plot_fwd_1_tomography_2D_002.png
:alt: Positions vs. Arrival Time
:srcset: /content/user-guide/tutorials/12-seismic/images/sphx_glr_plot_fwd_1_tomography_2D_002.png
:class: sphx-glr-single-img
.. rst-class:: sphx-glr-script-out
.. code-block:: none
<matplotlib.legend.Legend object at 0x7f9560aa9720>
.. GENERATED FROM PYTHON SOURCE LINES 162-167
Optional: Exporting Results
---------------------------
Write the data and true model
.. GENERATED FROM PYTHON SOURCE LINES 167-189
.. code-block:: Python
if save_file:
dir_path = os.path.dirname(tomo.__file__).split(os.path.sep)[:-3]
dir_path.extend(["tutorials", "seismic", "assets"])
dir_path = os.path.sep.join(dir_path) + os.path.sep
noise = 0.05 * dpred * np.random.randn(len(dpred))
data_array = np.c_[
np.kron(x, np.ones(n_receiver)),
np.kron(y_sources, np.ones(n_receiver)),
np.kron(np.ones(n_source), x),
np.kron(np.ones(n_source), y_receivers),
dpred + noise,
]
fname = dir_path + "tomography2D_data.obs"
np.savetxt(fname, data_array, fmt="%.4e")
output_model = model
fname = dir_path + "true_model_2D.txt"
np.savetxt(fname, output_model, fmt="%.4e")
.. rst-class:: sphx-glr-timing
**Total running time of the script:** (0 minutes 9.235 seconds)
**Estimated memory usage:** 288 MB
.. _sphx_glr_download_content_user-guide_tutorials_12-seismic_plot_fwd_1_tomography_2D.py:
.. only:: html
.. container:: sphx-glr-footer sphx-glr-footer-example
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download Jupyter notebook: plot_fwd_1_tomography_2D.ipynb <plot_fwd_1_tomography_2D.ipynb>`
.. container:: sphx-glr-download sphx-glr-download-python
:download:`Download Python source code: plot_fwd_1_tomography_2D.py <plot_fwd_1_tomography_2D.py>`
.. container:: sphx-glr-download sphx-glr-download-zip
:download:`Download zipped: plot_fwd_1_tomography_2D.zip <plot_fwd_1_tomography_2D.zip>`
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_