Commit a91128be authored by Tilman Steinweg's avatar Tilman Steinweg

updated manual

Still changes to be done in toy example chapter.
(figures with rotated coordinate system)
parent 330ad3e0
......@@ -10,9 +10,9 @@ The code was tested on our local workstation cluster using Suse Linux and OpenMP
\begin{small}
\begin{tabular}{lll}
Program & Description & Weblink \\
OpenMPI & MPI Implementation & \tiny{\url{http://www.mcs.anl.gov/research/projects/mpich2}} \\
OpenMPI & MPI Implementation & \url{http://www.open-mpi.org} \\
& (for the parallelization) & \\
MPICH2 & MPI Implementation & \url{http://www.open-mpi.org} \\
MPICH2 & MPI Implementation & \tiny{\url{http://www.mcs.anl.gov/research/projects/mpich2}} \\
& (for the parallelization) & \\
LAM/MPI & MPI Implementation & \url{http://www.lam-mpi.org} \\
& (for the parallelization) & \\
......@@ -37,7 +37,7 @@ After unpacking the software package (e.g. by \textbf{tar -zxvf ifos3D.tgz}) an
\textbf{bin}\\
This directory contains all executable programs, like ifos3D and snapmerge. These executables are generated using the command \textbf{make $<$program$>$} (see below).\vspace{0.2cm}\\
\textbf{doc}\\
This directory contains documentation on the software (this users guide) as well as some important papers in PDF format on which the software is based on.\vspace{0.2cm}\\
This directory contains documentation on the software (this users guide). \vspace{0.2cm}\\
\textbf{mfiles}\\
Here some Matlab routines (m-files) are stored. They can be used vor processing and visualisation of data. \vspace{0.2cm}\\
\textbf{par}\\
......@@ -59,10 +59,10 @@ This directory contains the complete source codes. The different subprograms are
\textbf{libcseife}\\
The libcseife library serves for filtering seismograms
\section{Compilation}
Before compiling the main program ifos3D, the libcseife library for the filtering of seismograms needs to be installed. In order to do this change to the directory libcseife. Remove the *.d files and perform \textbf{make all}. If problems with the compilation arise open \textit{Makefile} and adjust the compiler options for your system.\\
To compile the main program IFOS3D you can use the shell script \textit{compileIFOS3D.sh} in the \textit{par}-directory. This script executes make ifos3D in the \textit{src}-directory. To change the compiler options open the \textit{Makefile} in \textit{src}. Here you can also find examples for compiler options on different systems where IFOS3D or SOFI3D were used in the past.
To compile ifos3D change to the \textit{par}-directory and perform \textbf{make}. The Makefile will first compile the external libseife library and then compile the main program. In some cases it's necessary to remove the *.d files in the \textit{libseife}-directory. If problems with the compilation arise go to the \textit{libseife}-directory open \textit{Makefile} and adjust the compiler options for your system.\\
To change the compiler options open the \textit{Makefile} in \textit{src}-directory. Here you can also find examples for compiler options on different systems where IFOS3D or SOFI3D were used in the past.
\section{Running IFOS3D}
To start the program with OpenMPI you can use the command\\
\textbf{mpirun -np 8 nice -19 ../bin/ifos3D ./in\_and\_out/ifOS3D\_toy.inp $\mid$ tee ./in\_and\_out/ifos3D.out}\\
\textbf{mpirun -np 8 nice -19 ../bin/ifos3D ./in\_and\_out/ifOS3D\_toy.json $\mid$ tee ./in\_and\_out/ifos3D.out}\\
which runs IFOS3D on 8 processors with lowest priority. The standard output is written to \textit{/in\_and\_out/ifos3D.out}. You can also use the shell-script \textit{par/startIFOS3D.sh}.\\
For 3D FWI applications it is often useful to use supercomputers in high performance computing centers. Some examples for job scripts can be found in the directory \textit{scripts}.
......@@ -28,7 +28,8 @@ The elastic medium is described by the density $\rho$ and the Lam\'e parameters
\subsection{Finite difference modeling}\label{sec:FDmodeling}
For the simulations of 3D elastic wave propagation these equations are approximated by finite differences (FD). For the calculation of partial derivatives, the wavefield parameters and model parameters are discretised on a staggered-grid system \citep{Vir86, Lev88}. The 3D staggered-grid system is sketched in Figure~\ref{fig:stag_grid}. We use a Cartesian grid (x,y,z)=(i,j,k). The model parameters $\lambda$, $\mu$ and $\rho$ and the diagonal stress components $\sigma_{ii}$ are localised on the full grid points, whereas velocities and off-diagonal stress components are calculated on a grid which is half a grid point shifted to the original system. Compared to a conventional Cartesian grid, the staggered grid enables a larger grid spacing for the same level of accuracy.\\
\begin{figure}
\includegraphics[width=0.7\textwidth]{fig/staggered_grid}
\centering
\includegraphics[width=0.74\textwidth]{fig/ssg_grid}
\caption[Staggered grid system]{Staggered-grid coordinate system used for 3D FD modelling in SOFI3D}\label{fig:stag_grid}
\end{figure}
The derivatives are approximated by centered finite differences. The easiest finite-differences are of second order, where the two adjacent grid points with a grid space $dh$ are used to calculate the derivative between these grid points as follows
......@@ -46,7 +47,7 @@ Hereby the summation of the update values over all time steps accounts for the t
\subsubsection*{Boundary conditions}
There are two ways to avoid arteficial reflections from model boundaries in SOFI3D. The first exponentially damps waves within a boundary zone surrounding the model by multiplying the amplitudes with an exponentially decaying factor \citep{Cer85}. At least 30 gridpoints thickness at boundary zones is required. The second method is the implementation of convolutional perfectly matched layers (C-PMLs) \citep{Ber94,Kom07}. This method can offer a much better performance compared to the conventional exponential damping and a boundary zone of 10 gridpoints thickness is generally sufficient. However, in case of strong contrasts and heterogeneities at the model boundary instabilities can be caused.
\subsubsection*{The free surface}
At the free surface, the vertical stress components $\sigma_{xz}$, $\sigma_{yz}$ and $\sigma_{zz}$ are zero. A planar free surface is implemented implicitely in SOFI3D with the mirroring technique described by \cite{Lev88}.
At the free surface, the vertical stress components $\sigma_{xy}$, $\sigma_{yy}$ and $\sigma_{zy}$ are zero. A planar free surface is implemented implicitely in SOFI3D with the mirroring technique described by \cite{Lev88}.
\subsubsection*{Grid dispersion and stability}
To mitigate numerical dispersion and grid anisotropy of the wavefield the grid spacing must be chosen sufficiently small. For a fourth-order spatial and second-order temporal scheme the criterion
\begin{equation}
......
\chapter{In- and output}
\section{The input file - ifos3D.inp}
In the following we will describe the input parameters of IFOS3D. The grey text shows, how they appear in \textit{in\_and\_out/ifos3D.inp}. The parameters used for the wavefield simulation mostly resemble the parameters of SOFI3D and are also described in the SOFI manual. Note, that the format of the IFOS3D input file differs from SOFI3D.
\section{The input file - *.json}
In the following we will describe the input parameters of IFOS3D. The verbatim text shows, how they appear in \textit{in\_and\_out/ifos3d\_inv\_all\_parameters.json}. There is also a commented json file: \textit{in\_and\_out/ifos3d\_inv\_all\_parameters\_commented.json}. The parameters used for the wavefield simulation mostly resemble the parameters of SOFI3D and are also described in the SOFI manual.
\subsection{The grid and its decomposition}
\subsubsection*{The grid}
\textcolor{Gray}{number\_of\_gridpoints\_in\_x-direction\_(NX) = 160\\
number\_of\_gridpoints\_in\_y-direction\_(NY) = 160\\
number\_of\_gridpoints\_in\_z-direction\_(NZ) = 184\\
distance\_between\_gridpoints(in\_m)\_in\_x-direction\_(DX) = 0.8\\
distance\_between\_gridpoints(in\_m)\_in\_y-direction\_(DY) = 0.8\\
distance\_between\_gridpoints(in\_m)\_in\_z-direction\_(DZ) = 0.8}\vspace{0.1cm}\\
The cartesian grid system is specified by the number of grid points in each direction (\textbf{NX, NY, NZ}) and by the distance between grid points \textbf{DX, DY, DZ} chosen in the input file. Hereby \textbf{NZ} and \textbf{DZ} denote the vertical direction. This is different from the current SOFI3D release, where $Y$ denotes the vertical coordinate. Note, that internally the IFOS3D code also uses $y$ as vertical coordinate. \\
\begin{verbatim}
"Note that y denotes the vertical direction !" : "comment",
"3-D Grid" : "comment",
"NX" : "160",
"NY" : "184",
"NZ" : "160",
"DX" : "0.8",
"DY" : "0.8",
"DZ" : "0.8",
\end{verbatim}
The cartesian grid system is specified by the number of grid points in each direction (\textbf{NX, NY, NZ}) and by the distance between grid points \textbf{DX, DY, DZ} chosen in the input file. Hereby \textbf{NY} and \textbf{DY} denote the vertical direction!
The grid spacing needs to satisfy the dispersion criterion at least up to the maximum frequency and minimum wavelength used in the inversion. In case the dispersion criterion is violated for the frequency range of the source wavelet, a warning is given at the beginning of the simulation. The dispersion criterions are given in the SOFI3D manual.
\subsubsection*{Domain decomposition}
\textcolor {Gray}{number\_of\_processors\_in\_x-direction\_(NPROCX) = 8\\
number\_of\_processors\_in\_y-direction\_(NPROCY) = 8\\
number\_of\_processors\_in\_z-direction\_(NPROCZ) = 8} \vspace{0.1cm}\\
\begin{verbatim}
"Domain Decomposition" : "comment",
"NPROCX" : "2",
"NPROCY" : "2",
"NPROCZ" : "2",
\end{verbatim}
The parallelisation is based on domain decomposition. Each processor performs the FWI algorithm for a small subvolume of the total grid system. \textbf{NPROCX}, \textbf{NPROCY} and \textbf{NPROCZ} define the number of processors in each direction, which gives a total number of $NPROCX\times NPROCY\times NPROCZ$ processors. Note, that the number of processors in each direction needs to be a common factor of the number of grid points in this direction.
\subsection{FD-modeling parameters}
\subsubsection*{Order of FD operator}
\textcolor{Gray}{
order\_of\_spatial\_fd\_operators\_(FDORDER) = 4\\
fd\_coefficients\_(Taylor=1;Holberg=2)\_(FDCOEFF) = 2}\vspace{0.1cm}\\
The order of the spatial FD-sheme, which is used for the stress and velocity updates can be chosen as 2, 4, 6, 8, 10, and 12. These orders result in different dispersion and stability criterions (section~\ref{sec:SOFI3Dcharacter}). A larger FD-order enables a larger grid spacing, however in case of strong heterogeneities a smaller FD-operatot is more accurate. Within the CPML boundary a fourth order operator is applied automatically. We also use fourth order operators to calculate the spatial derivatives of the wavefield velocity, which are needed for the gradient calculation. IFOS3D uses a second order time discretisation.
\begin{verbatim}
"FD order" : "comment",
"FDORDER" : "4",
"FDCOEFF" : "2",
\end{verbatim}
The order of the spatial FD-sheme, which is used for the stress and velocity updates can be chosen as \textbf{FDORDER}=2, 4, 6, 8, 10, and 12.
These orders result in different dispersion and stability criterions (section~\ref{sec:SOFI3Dcharacter}). A larger FD-order enables a larger grid spacing, however in case of strong heterogeneities a smaller FD-operatot is more accurate.
With the option \textbf{FDCOEFF} the user can switch between Taylor (\textbf{FDCOEFF}=1) and Holberg (\textbf{FDCOEFF}=2) FD coefficients. Within the CPML boundary a fourth order operator is applied automatically.
We also use fourth order operators to calculate the spatial derivatives of the wavefield velocity, which are needed for the gradient calculation. IFOS3D uses a second order time discretisation.
\subsubsection*{Time Stepping}
\textcolor {Gray}{
time\_of\_wave\_propagation\_(in\_sec)\_(TIME) = 0.06\\
timestep\_(in\_seconds)\_(DT) = 5.0}\vspace{0.1cm}\\
TIME defines the total simulation length. The time stepping DT must be chosen to satisfy the stability criterion (SOFI3D maunual) and thus depends on the grid spacing and the maximum velocity. In case of violation IFOS3D will give a warning and terminate. The number of timesteps of the inversion is given by $TIME/DT$.
\begin{verbatim}
"Time Stepping" : "comment",
"TIME" : "0.06",
"DT" : "5.0e-05",
\end{verbatim}
\textbf{TIME} defines the total simulation length. The time stepping \textbf{DT} must be chosen to satisfy the stability criterion (SOFI3D maunual) and thus depends on the grid spacing and the maximum velocity.
In case of violation IFOS3D will give a warning and terminate. The number of timesteps of the inversion is given by $TIME/DT$.
\subsection{Sources}
\begin{verbatim}
"Source" : "comment",
"SOURCE_SHAPE" : "4",
"SOURCE_TYPE" : "4",
"ALPHA" : "45.0",
"BETA" : "45.0",
"SIGNAL_FILE" : "./STF/stf.su",
"SRCREC" : "1",
"SOURCE_FILE" : "./sources/sources_toy.dat",
"RUN_MULTIPLE_SHOTS" : "1",
\end{verbatim}
\subsubsection*{The wavelet}
\textcolor {Gray}{(ricker=1;fumue=2;from\_SIGNAL\_FILE=3;SIN**3=4)\_(QUELLART) = 4\\
SIGNAL\_FILE = ?} \vspace{0.1cm}\\
Different wavelets for the modeling can be chosen:
\begin{itemize}
\item Ricker wavelet: $s=(1.0-2.0\tau ^2)\text{exp}(-\tau ^2); \hspace{0.5cm} \tau=\pi f_c(t-1.5/f_c-t_d);$
\item Fuchs-M\"uller wavelet: $s(t)=\text{sin}(2.0\pi (t-t_d)f_c-0.5 \text{sin}(4.0\pi (t-t_d)f_c); \hspace{0.2cm}\text{for}\hspace{0.2cm}(t_d<t<t_d+1/f_c) \hspace{0.5cm}\text{else}\hspace{0.2cm} s(t)=0$
\item sin$^3$-wavelet: $s(t)=0.75\pi f_c\text{sin}(\pi f_c(t-t_d)^3;\hspace{0.2cm}\text{for}\hspace{0.2cm}(t_d<t<t_d+1/f_c) \hspace{0.5cm}\text{else}\hspace{0.2cm} s(t)=0$
\item step function: $s(t)=0.0; \hspace{0.2cm}\text{if}\hspace{0.2cm}(t<t_d) \hspace{0.5cm}\text{else} \hspace{0.2cm} s(t)=1.0; $
\item \textbf{SOURCE\_SHAPE}=1: Ricker wavelet:\\ $s=(1.0-2.0\tau ^2)\text{exp}(-\tau ^2); \hspace{0.5cm} \tau=\pi f_c(t-1.5/f_c-t_d);$
\item \textbf{SOURCE\_SHAPE}=2: Fuchs-M\"uller wavelet:\\ $s(t)=\text{sin}(2.0\pi (t-t_d)f_c-0.5 \text{sin}(4.0\pi (t-t_d)f_c); \hspace{0.2cm}\text{for}\hspace{0.2cm}(t_d<t<t_d+1/f_c) \hspace{0.4cm}\text{else}\hspace{0.2cm} s(t)=0$
\item \textbf{SOURCE\_SHAPE}=4: sin$^3$-wavelet:\\ $s(t)=0.75\pi f_c\text{sin}(\pi f_c(t-t_d)^3;\hspace{0.2cm}\text{for}\hspace{0.2cm}(t_d<t<t_d+1/f_c) \hspace{0.5cm}\text{else}\hspace{0.2cm} s(t)=0$
\item \textbf{SOURCE\_SHAPE}=5: step function:\\ $s(t)=0.0; \hspace{0.2cm}\text{if}\hspace{0.2cm}(t<t_d) \hspace{0.5cm}\text{else} \hspace{0.2cm} s(t)=1.0; $
\end{itemize}
The employed center frequency $f_c$ and the time delay of the source wavelet $t_d$ are defined in the source file. An exemplary plot of the different wavelets and the corresponding spectra can be seen in the SOFI3D manual. Note, that the symmetric Ricker wavelet is delayed and its maximum period is excited at the source location after one period. The step function is used for the calculation of the diagonal Hessian approximation. \\
If you want to define your own wavelet, you can include this either in the source code in \textit{src/wavelet.c} or use QUELLART=3 and read in an external wavelet from SIGNAL\_FILE. This file should contain the source signal in ASCII with the correct time sampling and with one sample per line, like \\
The employed center frequency $f_c$ and the time delay of the source wavelet $t_d$ are defined in the source file. An exemplary plot of the different wavelets and the corresponding spectra can be seen in the SOFI3D manual.
Note, that the symmetric Ricker wavelet is delayed and its maximum period is excited at the source location after one period. The step function is used for the calculation of the diagonal Hessian approximation. \\
If you want to define your own wavelet, you can include this either in the source code in \textit{src/wavelet.c} or use \textbf{SOURCE\_SHAPE}=3 and read in an external wavelet from \textbf{SIGNAL\_FILE}.
This file should contain the source signal in ASCII with the correct time sampling and with one sample per line, like \\
----------------\\
0.0\\
0.01\\
......@@ -49,55 +80,73 @@ If you want to define your own wavelet, you can include this either in the sourc
...\\
----------------
\subsubsection*{The source type}
\textcolor {Gray}{point\_source\_(explosive=1;force\_in\_x=2;in\_y=3;in\_z=4;custom=5)\_(QUELLTYP) = 4\\
force\_angle\_between\_x\_y\_(in\_degree)\_(ALPHA) = 45.0\\
force\_angle\_between\_x\_z\_(in\_degree)\_(BETA) = 45.0}\vspace{0.1cm}\\
You can either choose an explosive source exciting compressional waves or a point force directed in $x$, $y$ or $z$ (vertical) direction. With QUELLTYP=5 it is also possible to define a force in arbitrary direction by using the angles ALPHA and BETA. This is illustrated in the SOFI manual. The plane wave excitation was not yet tested for inversion in IFOS3D but can also be chosen here.
You can either choose an explosive source (\textbf{SOURCE\_TYPE}=1) exciting compressional waves or a point force directed in $x$, $y$ (vertical) or $z$ (\textbf{SOURCE\_TYPE}=2,3,4) direction. With \textbf{SOURCE\_TYPE}=5 it is also possible to define a force in arbitrary direction by using the angles \textbf{ALPHA} and \textbf{BETA}. This is illustrated in the SOFI manual.
The plane wave excitation (\textbf{SOURCE\_REC}=2) was not yet tested for inversion in IFOS3D but can also be chosen here.
%# Plane wave excitation,if PLANE_WAVE_DEPTH>0, SRCREC is treated as 0
%depth_of_plane_wave_excitation_(no<=0)_(in_meter)_(PLANE_WAVE_DEPTH) = 0.0
%dip_of_plane_wave_from_vertical_(in_degrees)_(PHI) =0.0
%duration_of_source-signal_(in_seconds)_(TS) = 0.0033
\subsubsection*{The source location}
\textcolor {Gray}{read\_source\_positions\_from\_SOURCE\_FILE\_(yes=1)\_(SRCREC) = 1\\
SOURCE\_FILE = ./sources/sources\_florian.dat \\
run\_multiple\_shots\_defined\_in\_SOURCE\_FILE\_(yes=1)\_(RUN\_MULTIPLE\_SHOTS) = 1}\vspace{0.1cm}\\
The source locations are defined in the SOURCE\_FILE located in the folder \textit{sources}. For each source the parameters XSRC, YSRC and ZSRC are defined, which are the $x$-, $y$- and $z$-coordinates of the source position in meter. Additionally a time delay TD in seconds, a center frequency of the source signal (FC) in Hz and an amplitude of the source signal are chosen. An example file for a source located at (XSRC=100m, YSRC=50m and ZSRC=2m) with a center frequency of 20Hz and no time delay looks like:\\
The source locations are defined in the \textbf{SOURCE\_FILE} located in the folder \textit{sources}. For each source the parameters XSRC, YSRC and ZSRC are defined, which are the $x$-, $y$- and $z$-coordinates of the source position in meter. Additionally a time delay TD in seconds, a center frequency of the source signal (FC) in Hz and an amplitude of the source signal are chosen. An example file for a source located at (XSRC=100m, YSRC=50m and ZSRC=2m) with a center frequency of 20Hz and no time delay looks like:\\
----------------\\
100.0 \hspace{0.2cm} 50.0 \hspace{0.2cm} 2.0 \hspace{0.2cm} 0.0 \hspace{0.2cm} 20.0 \hspace{0.2cm} 1.0\\
----------------
\subsection{The model}
\subsubsection*{Model input}
\textcolor {Gray}{
read\_model\_from\_MFILE(yes=1)(READMOD) = 0\\
MFILE = model/toy} \vspace{0.1cm}\\
\begin{verbatim}
"Model" : "comment",
"READMOD" : "0",
"MFILE" : "model/toy",
\end{verbatim}
The elastic subsurface model is parametrised by the compressional wave velocity $v_p$, the shear wave velocity $v_s$ and the density $\rho$ at each grid point. Model parameters are required for the forward modeling and as starting model for the inversion. There are two options to read in the model parameters:\\
READMOD=1: The model parameters of an external model are read in from the model files defined as MFILE, here for example \textit{toy.$v_p$}, \textit{toy.$v_s$} and \textit{toy.$\text{rho}$} placed in the folder \textit{par/model}. In these files, each material parameter value must be saved as 32 bit (4 byte) native float.
\textbf{READMOD}=1: The model parameters of an external model are read in from the model files defined as \textbf{MFILE}, here for example \textit{toy.$v_p$}, \textit{toy.$v_s$} and \textit{toy.$\text{rho}$} placed in the folder \textit{par/model}. In these files, each material parameter value must be saved as 32 bit (4 byte) native float.
Velocities must be in m/s, density values in kg/m$^3$ defined at each grid point. The values are read in with \textit{src/readmod.c} and you may look at this source code for the correct order of parameters.\\
READMOD=0: It is also possible to generate the model parameters in the program on the fly. The corresponding C-function, e.g. \textit{src/hh.c} is included in the \textit{src/Makefile}.
\textbf{READMOD}=0: It is also possible to generate the model parameters in the program on the fly. The corresponding C-function, e.g. \textit{src/hh.c} is included in the \textit{src/Makefile}.
\subsubsection*{Viscoelasticity}
\textcolor {Gray}{Number\_of\_relaxation\_mechanisms\_(L) = 0\\
L\_Relaxation\_frequencies\_(FL) = 1000.0\\
Tau\_value\_for\_entire\_model\_(TAU) = 0.000001}\vspace{0.1cm}\\
The inversion code was only tested in the elastic version (L=0). However, the viscoelasic version of SOFI3D is included in the package and can be used for forward modeling. In this case, the elastic update functions for the velocities need to be replaced by the viscoelastic update functions available in \textit{src}. No inversion code is included for viscoelastic parameters and they can be only used as passive parameters. For further information on viscoelasticity please look in the SOFI manual.
\begin{verbatim}
"Q-approximation" : "comment",
"L" : "0",
"FL1" : "1000.0",
"TAU" : "0.000001",
\end{verbatim}
The inversion code was only tested in the elastic version (\textbf{L}=0). However, the viscoelasic version of SOFI3D is included in the package and can be used for forward modeling. In this case, the elastic update functions for the velocities need to be replaced by the viscoelastic update functions available in \textit{src}. No inversion code is included for viscoelastic parameters and they can be only used as passive parameters. For further information on viscoelasticity please look in the SOFI manual.
\subsection{Boundary conditions}
\subsubsection*{The free surface}
\textcolor {Gray}{free\_surface\_(yes=1)(FREE\_SURF) = 0}\vspace{0.1cm}\\
\begin{verbatim}
"Free Surface" : "comment",
"FREE_SURF" : "0",
\end{verbatim}
A plane stress free surface is applied at the top of the global grid if FREE SURF = 1 using the imaging method proposed by \citep{Lev88}. If FREE SURF = 0 a full space is simulated.
\subsubsection*{Model boundaries}
\textcolor {Gray}{type\_of\_boundary\_condition\_(ABS\_TYPE)\_(PML=1/ABS=2) = 1\\
width\_of\_absorbing\_frame\_(in\_grid\_points)\_(No$<$=0)\_(FW) = 10\\
percentage\_of\_amplitude\_decay\_at\_outer\_edge\_(DAMPING) = 8.0\\
Dominant\_frequency\_(in\_Hz)\_(FPML) = 200.0\\
P\_wave\_velocity\_near\_the\_grid\_boundary(in\_m/s)\_(VPPML) = 6200.0\\
apply\_periodic\_boundary\_condition\_at\_edges\_(BOUNDARY):\\
(no=0)\_(left/right/front/back=1) = 0 }\vspace{0.1cm}\\
\begin{verbatim}
"Absorbing Boundary" : "comment",
"ABS_TYPE" : "1",
"FW" : "10",
"DAMPING" : "8.0",
"VPPML" : "6200.0",
"FPML" : "200.00",
"BOUNDARY" : "0",
\end{verbatim}
The model includes a boundary layer of \textbf{FW} gridpoints at each side, where the wavefield is damped. There are two options included in IFOS3D:\\
\textbf{ABS\_TYPE}=1: IFOS3D offers the use of convolutional perfectly matching (C-PML) layers to damp wavefields in the model boundaries. This technique shows a very good ability to remove reflections from the model boundaries. The parameters \textbf{FPML} as the dominant frequency and \textbf{VPPML} as the compressional wave velocity near the boundary are employed to calculate the C-PML coefficients. For the C-PML techique a width of \textbf{FW}=10 gridpoints is generally enough. Sometimes, instabilities can be created in the PML boundaries, especially in case of strong model heterogeneities in this area. For more information we refer to the SOFI manual.\\
\textbf{ABS\_TYPE}=2: The ``traditional`` way to prevent reflections from the model boundary is the exponential damping. This technique is very robust, but less efficient. A thickness of at least \textbf{FW}=30 gridpoints is required. For the parameter \textbf{DAMPING} you can use about 8\%. Note that much higher values cause reflections at the boundary interface.\\
The periodic boundary condition (\textbf{BOUNDARY}=1) is taken from SOFI3D. It was not yet tested in the IFOS3D code and needs to be fully implemented and tested before use.
\subsection{Seismogram output}
\textcolor {Gray}{output\_of\_seismograms\_(SEISMO) = 1}\vspace{0.1cm}\\
\begin{verbatim}
"Receiver" : "comment",
"SEISMO" : "1",
"READREC" : "0",
"REC_FILE" : "./receiver/receiver.dat",
"REFRECX, REFRECY, REFRECZ" : "0.0 , 0.0, 0.0",
"XREC1, YREC1, ZREC1" : "90.0 , 90.0, 90.0",
"XREC2, YREC2, ZREC2" : "90.0 , 90.0, 90.0",
"NGEOPH" : "1",
\end{verbatim}
In every iteration and for each shot IFOS3D stores the seimograms. The parameter \textbf{SEISMO} defines the output:
\begin{itemize}
\item 0: no seismograms
......@@ -106,30 +155,41 @@ In every iteration and for each shot IFOS3D stores the seimograms. The parameter
\item 3: curl and div
\item 4: everything
\end{itemize}
For the inversion with IFOS3D we generally output particle velocities (one file for each component), which can be used to look at the fit of the waveforms. The inversion of pressure wavefields was not tested with IFOS3D. For more information about the curl and div output we refer to the SOFI manual.
For the inversion with IFOS3D we generally output particle velocities (one file for each component), which can be used to look at the fit of the waveforms. The inversion of pressure wavefields was not tested with IFOS3D.
For more information about the curl and div output we refer to the SOFI manual.
\subsubsection*{Receiver locations}
The receivers define the location of the seismogram output. They are always located on the grid points (no interpolation) and if necessary they are shifted to the nearest grid point. Note, that $z$ defines the vertical direction. There are three options to define receiver locations in IFOS3D:\vspace{0.3cm}\\
\textcolor {Gray}{read\_receiver\_positions\_from\_file\_(yes=1)\_(READREC) = 0\\
REC\_FILE = ./receiver/receiver.dat\\
reference\_point\_for\_receiver\_coordinate\_system\_(REFREC) = 0.0 , 0.0 , 0.0}\vspace{0.1cm}\\
If the option \textbf{READREC}=1 is used, receivers are read from \textbf{REC\_FILE}. This ASCII file defines one receiver in a row using the $x-$, $y$- and $z$ coordinates in meter. It is possible to shift the receiver locations by employing a reference point defined by \textbf{REFREC}.
\vspace{0.3cm}\\
\textcolor {Gray}{position\_of\_first\_receiver\_(in\_m)\_(XREC1,YREC1,ZREC1) = 90.0 , 90.0 , 90.0\\
position\_of\_last\_receiver\_(in\_m)\_(XREC2,YREC2,ZREC2) = 90.0 , 90.0 , 90.0\\
distance\_between\_two\_adjacent\_receivers\_(in\_gridpoints)\_(NGEOPH) = 1}\vspace{0.1cm}\\
The second option (\textbf{READREC}=0) defines a straight line of receivers, like a horizontal or vertical profile. This profile is defined by its starting point (\textbf{XREC1,YREC1,ZREC1}), its end point (\textbf{XREC2,YREC2,ZREC2}) and the distance between the receivers (\textbf{NGEOPH}).
\vspace{0.3cm}\\
\textcolor {Gray}{number\_of\_planes\_(no$<$=0)\_(REC\_ARRAY) = 1\\
depth\_of\_first\_(upper)\_plane\_(in\_m)\_(REC\_ARRAY\_DEPTH) = 24.0\\
vertical\_distance\_between\_planes\_(in\_m)\_(REC\_ARRAY\_DIST) = 30.0 \\
distance\_between\_receivers\_in\_x-direction\_(in\_gridpoints)\_(DRX) = 10\\
distance\_between\_receivers\_in\_y-direction\_(in\_gridpoints)\_(DRY) = 10}\vspace{0.1cm}\\
The last option is very attractive for 3D FWI because it defines seismic arrays which are very useful for 3D synthetic FWI applications. It is possible to define several horizontal arrays (\textbf{REC\_ARRAY}) in different depths. The upper plane is at \textbf{REC\_ARRAY\_DEPTH} with a vertical distance of \textbf{REC\_ARRAY\_DIST} in meter to the next receiver plane. The horizontal distance between receivers is defined by \textbf{DRX} and \textbf{DRY}. Note, that the receiver array starts in a distance of 10 gridpoints from the absorbing boundary.
The receivers define the location of the seismogram output. They are always located on the grid points (no interpolation) and if necessary they are shifted to the nearest grid point. Note, that $y$ defines the vertical direction.
There are three options to define receiver locations in IFOS3D:\vspace{0.3cm}\\
If the option \textbf{READREC}=1 is used, receivers are read from \textbf{REC\_FILE}. This ASCII file defines one receiver in a row using the $x-$, $y$- and $z$ coordinates in meter.
It is possible to shift the receiver locations by employing a reference point defined by \textbf{REFREC}.
The second option (\textbf{READREC}=0) defines a straight line of receivers, like a horizontal or vertical profile. This profile is defined by its starting point (\textbf{XREC1,YREC1,ZREC1}),
its end point (\textbf{XREC2,YREC2,ZREC2}) and the distance between the receivers (\textbf{NGEOPH}).
\begin{verbatim}
"Receiver array" : "comment",
"REC_ARRAY" : "1",
"REC_ARRAY_DEPTH" : "24.0",
"REC_ARRAY_DIST" : "30.0",
"DRX" : "10",
"DRZ" : "10",
\end{verbatim}
The last option (\textbf{READREC}=2) is very attractive for 3D FWI because it defines seismic arrays which are very useful for 3D synthetic FWI applications. It is possible to define several horizontal arrays (\textbf{REC\_ARRAY}) in different depths. The upper plane is at \textbf{REC\_ARRAY\_DEPTH} with a vertical distance of \textbf{REC\_ARRAY\_DIST} in meter to the next receiver plane. The horizontal distance between receivers is defined by \textbf{DRX} and \textbf{DRY}. Note, that the receiver array starts in a distance of 10 gridpoints from the absorbing boundary.
\subsubsection*{Seismograms}
\textcolor {Gray}{samplingrate\_and\_timelag\_(in\_timesteps!)\_(NDT,NDTSHIFT) = 1 , 0\\
data-format\_(SEIS\_FORMAT[6]) = 1\\
basic\_filename\_(SEIS\_FILE) = ./su/cal\_toy}\vspace{0.1cm}\\
\begin{verbatim}
"Seismograms" : "comment",
"NDT" : "1",
"NDTSHIFT" : "0",
"SEIS_FORMAT" : "1",
"SEIS_FILE" : "su/IFOS",
\end{verbatim}
For the seismogram output every \textbf{NDT}th sample is written to file starting at timestep \textbf{NDTSHIFT} of the FD modeling. The file format can be choosen with \textbf{SEIS\_FORMAT}:\\
0: SEG-Y (ASCII-text/native 4-byte-floats (IEEE on PC)/little endian on PC)\\
1: SU (native 4-byte-floats (IEEE on PC)/little endian on PC)\\
......@@ -140,33 +200,53 @@ For the seismogram output every \textbf{NDT}th sample is written to file startin
We recommend the use of the SU-format. This format stores a header in front of each trace with information like time sampling, source and receiver position and tracenumber. Seismic Unix offers different functions for data processing of SU-files.\\
The seismograms are stored for each component, each shot and each iteration in the folder \textit{su}. \textbf{SEIS\_FILE} defines the name of the files, here for example \textit{cal\_toy\_vx\_it3.su.shot2}
\subsubsection*{Snapshots}
\textcolor {Gray}{output\_of\_snapshots\_(SNAP)(yes$>$0) = 0}\vspace{0.1cm}\\
\begin{verbatim}
"Snapshots" : "comment",
"SNAP" : "0",
\end{verbatim}
It is also possible to store snapshots of the wavefield (\textbf{SNAP}$>$0), that means the wavefield in the whole volume for a series of times. This output is generally not useful in a full inversion due to the great amount of stored data. It can however be used for the forward modeling. It is also possible to save snapshots of the backpropagated wavefield, however in this case the function \textit{snap()} needs to be included in the timeloop of the backpropagation. For the different options to save snapshots and for the processing (merging and visualisation) of the snapshots we refer to the SOFI3D manual.
\subsection{Forward modeling and inversion}
\textcolor {Gray}{method\_(METHOD) = 1}\vspace{0.1cm}\\
\begin{verbatim}
"Method" : "comment",
"METHOD" : "0",
\end{verbatim}
Using \textbf{METHOD}=0 IFOS3D performs a forward simulation similar to SOFI3D. For synthetic tests this option can be used to create observed data of the ''true`` model. For this option the stored seismograms are stored for each processor containing a receiver seperately and as a merged file for all processors. This way the data can be directly used as observed data in the FWI.\\
With the choice of \textbf{METHOD}=1 IFOS3D performs a conjugate gradient FWI as specified by the following parameters in the input file.
\subsection{General inversion parameters}
\textcolor {Gray}{minimum/maximum\_iteration\_number\_($>$0)\_(ITMIN,ITMAX) = 1 , 80 }\vspace{0.1cm}\\
\begin{verbatim}
"General" : "comment",
"ITMIN, ITMAX" : "1 , 80",
"FILT" : "1",
"NFMAX" : "5",
"TAST" : "100",
"VP0, VS0, RHO0" : "6200.0, 3600.0, 2800.0",
"WEIGHT_VP,WEIGHT_VS,WEIGHT_RHO" : "1.0, 1.0, 0.0",
\end{verbatim}
The parameters \textbf{ITMIN} and \textbf{ITMAX} define the total number of iterations performed by IFOS3D. Each inversion starts with \textbf{ITMIN}=1. However it is possible to split the inversion in parts and use \textbf{ITMIN, ITMAX}=(1,15) in a first step and then proceed by choosing \textbf{ITMIN, ITMAX}=(16,30).
\vspace{0.3cm}\\
\textcolor {Gray}{filtering\_(yes=1/no=0)\_(FILT) = 1}\vspace{0.1cm}\\
The gradients are only calculated for discrete frequencies, which works as a natural filter for the inversion with IFOS3D. However, the misfit is calculated in time domain to include the full waveforms. Thus it is reasonable to use filtered seismograms. By choosing \textbf{FILT}=1 IFOS3D filters the source wavelet and the observed seismograms below the maximum frequency of the current iteration stage with a Butterworth filter of 4th order. This way the misfit is calculated for data up to the maximum inversion frequency.
\vspace{0.3cm}\\
\textcolor {Gray}{maximum\_number\_frequencies\_per\_iteration\_(NFMAX) = 5}\vspace{0.1cm}\\
The parameter \textbf{NFMAX} defines the maximum number of frequencies employed in one inversion stage.
\vspace{0.3cm}\\
\textcolor {Gray}{number\_of\_timesteps\_per\_period\_used\_for\_inversion\_(TAST) = 100}\vspace{0.1cm}\\
The parameter (\textbf{TAST}) defines the time sampling used to extract the monochromatic frequency wavefields via discrete Fourier transformation. The period used for the calculation of the sample rate corresponds to $1/f_{max}$, where $f_{max}$ is the maximum frequency of the current frequency stage. In general it is not necessary to use the small FD time sampling for the discrete Fourier transformation and runtime can be saved. However, when the number of timesteps used for the discrete FD-transformation is too low, the gradient looks pixelated. It is therefore recommendable to check the gradients when using this option.
\vspace{0.3cm}\\
\textcolor {Gray}{average\_model\_parameter\_(VP0,VS0,RHO0) = 6200.0, 3600.0, 2800.0}\vspace{0.1cm}\\
The parameters \textbf{VP0,VS0} and \textbf{RHO0} are reference values used for the different parameter classes, like for example the average of the starting model. They are used for the model update of the conjugate gradient method and for the parameter normalisation in the L-BFGS method.
\vspace{0.3cm}\\
\textcolor {Gray}{parameter\_class\_weighting\_factors\_for\_vp/vs/rho\_(WEIGHT) = 1.0, 1.0, 0.0}\vspace{0.1cm}\\
It is possible to employ weighting factors (\textbf{weight}) for each parameter class. A value of 1.0 means, that the parameter is fully updated. A value of 0.0 corresponds to no update of this parameter. This enables an empirical weighting of the different parameter classes.
It is possible to employ weighting factors (\textbf{WEIGHT}) for each parameter class. A value of 1.0 means, that the parameter is fully updated. A value of 0.0 corresponds to no update of this parameter. This enables an empirical weighting of the different parameter classes.
\subsection{Observed data input}
\textcolor {Gray}{observed\_data\_fileneame\_(SEIS\_OBS\_FILE) = ./su\_obs/obs\_toy\\
external\_or\_internal\_observed\_data\_(EXTOBS) = 0}\vspace{0.1cm}\\
\begin{verbatim}
"SEIS_OBS_FILE" : "./su_obs/obs_toy",
"EXTOBS" : "0",
\end{verbatim}
The inversion requires observed or measured data as input, which is read in from the folder \textit{su\_obs}. These data are required to
\begin{itemize}
\item be in SU format
......@@ -178,11 +258,14 @@ For \textbf{EXTOBS}=0 the observed data is data produced by a forward modeling o
For \textbf{EXTOBS}=1 the observed data is read in from one file for all processors. In this case IFOS3D splits the data before starting the inversion and saves it to seperate files for each processor. When using this option the order of the seismograms is required to resemble the order of receivers in the receiver file.
\subsection{Gradients}
\subsubsection*{Gradient output}
\textcolor {Gray}{gradient\_filename\_(GRAD\_FILE) = ./grad/toy\_grad}\vspace{0.1cm}\\
\begin{verbatim}
"GRAD_FILE" : "./grad/toy_grad",
"DAMPTYPE" : "2",
\end{verbatim}
Gradients of the L2-norm based data misfit are calculated for each parameter class ($v_p$, $v_s$ and $\rho$) for each shot and summed over all shots and used frequencies. In each iteration raw gradients are stored for each grid point in the folder \textit{grad} in binary format. The order of the values resembles the order in the model files. Their file name include the label given in the input file (GRAD\_FILE), the maximum inversion frequency of this iteration, the parameter class and the iteration number (e.g. \textit{grad/toy\_grad.vp\_200.00Hz\_it6}). Additionally, the conjugate gradients are stored. They are labelled with numbers iteration+1000 (e.g.\textit{grad/toy\_grad.vp\_200.00Hz\_it1006}). \\
It is possible to save additional gradients by including the function \textit{savegrad()} at the required position in the \textit{src/ifos3d.c} source code. This can be useful to look at gradients at single shots or to test a preconditioning function.
\subsubsection*{Gradient preconditioning}
\textcolor {Gray}{Type\_of\_preconditioning\_(DAMPTYPE) = 2}\vspace{0.1cm}\\
The option \textbf{DAMPTYPE} defines a local preconditioning of the gradients at source and receiver positions as described in section~\ref{sec:preconditioning}. The following options are avalable:\\
0: no damping\\
1: Gaussian taper around receivers\\
......@@ -190,15 +273,25 @@ The option \textbf{DAMPTYPE} defines a local preconditioning of the gradients at
3: Gaussian tapering of source and receiver planes\\
Note, that additionally a taper of the PML boundaries is applied. If the preconditioning is not sufficient for your gradient, it is possible to change the taper parameters in \textit{src/precongrad.c}.
\subsection{Step length calculation}
\textcolor {Gray}{number\_of\_shots\_used\_for\_steplength\_estimation\_(NSHOTS\_STEP) = 4\\
initial\_test\_steplength\_(TESTSTEP) = 0.02}\vspace{0.1cm}\\
\begin{verbatim}
"Steplength estimation" : "comment",
"NSHOTS_STEP" : "4",
"TESTSTEP" : "0.02",
\end{verbatim}
The step length calculation is performed using the parabola method as decribed in section~\ref{sec:steplength}. To avoid long runtimes for the calcultion only a subset of \textbf{NSHOTS\_STEP} is used for the estimation of the optimal steplength. The parameter \textbf{TESTSTEP} defines an initial test steplength used in the first iteration of each frequency stage. In later iterations this test steplength is adjusted to the inversion proceedings. For a conjugate gradient inversion a value of about 0.02 is generally a good choice. This means that for steplength estimation an update of 2\% and 4\% of the average model parameter is tested and that a maximum model update of five percent is possible.
\subsection{Model output}
\textcolor {Gray}{model\_output\_filename\_(MOD\_OUT\_FILE) = ./model/toy}\vspace{0.1cm}\\
\begin{verbatim}
"MOD_OUT_FILE" : "./model/toy",
\end{verbatim}
After model update, the model parameters at each grid point are written to disc into the folder \textit{model}. One file for each parameter ($v_p$, $v_s$ and $\rho$) is written in each iteration. The structure of the binary files is similar to the gradient output. The filename is labelled by \textbf{MOD\_OUT\_FILE} (e.g. \textit{model/toy.vs\_it6}).
\subsection{The workflow}
\textcolor {Gray}{inversion\_parameter\_file\_(INV\_FILE) = ./in\_and\_out/workflow\_toy.dat}\vspace{0.1cm}\\
The workflow file (\textit{in\_and\_out/workflow.dat}) organises the different inversion stages. An example is given in figure~\ref{fig:workflow_input}. One frequency stage is defined in each row. This includes the number of iterations per stage, the number of frequencies and a list of these frequencies. We did not yet include an abortion criteria, so that the number of iterations in one stage is a fixed value and corresponds to the maximum numver of iterations given in the workflow. Note that the maximum number of frequencies used in one stage is defined as \textbf{NFMAX} in \textit{ifos3D.inp}. IFOS3D always starts to read the upper row. When an inversion is split up, the upper row must resemble the current frequency stage.\\
\begin{verbatim}
"INV_FILE" : "./in_and_out/workflow_toy.dat",
\end{verbatim}
The workflow file (\textit{in\_and\_out/workflow.dat}) organises the different inversion stages. An example is given in figure~\ref{fig:workflow_input}. One frequency stage is defined in each row. This includes the number of iterations per stage, the number of frequencies and a list of these frequencies. We did not yet include an abortion criteria, so that the number of iterations in one stage is a fixed value and corresponds to the maximum numver of iterations given in the workflow. Note that the maximum number of frequencies used in one stage is defined as \textbf{NFMAX} in \textit{*.json}. IFOS3D always starts to read the upper row. When an inversion is split up, the upper row must resemble the current frequency stage.\\
There are some parameters included in the workflow which are useful in FWI, but not yet implemented in IFOS3D. These include parameters for time and offset windowing and an abortion criteria. They will be included in future work.
\begin{figure}[h!]
\begin{center}
......@@ -207,20 +300,28 @@ There are some parameters included in the workflow which are useful in FWI, but
\end{center}
\end{figure}
\subsection{Hessian}
\textcolor {Gray}{Apply\_Hessian\_(yes=1/no=0)\_(HESS) = 0\\
Read\_Hessian\_from\_file\_(READ\_HESS) = 0\\
Part\_of\_receivers\_used\_for\_Hessian\_(REC\_HESS) = 1\\
Water\_level\_Hessian\_for\_vp/vs/rho\_(WATER\_HESS) = 0.0192, 0.0192, 1.0e-14\\
hessian\_file\_(HESS\_FILE) = ./hess/toy\_hess }\vspace{0.1cm}\\
\begin{verbatim}
"Hessian" : "comment",
"HESS" : "0",
"READ_HESS" : "0",
"REC_HESS" : "1",
"WATER_HESS_VP, WATER_HESS_VS, WATER_HESS_RHO" : "0.0192, 0.0192, 1.0e-14",
"HESS_FILE" : "./hess/toy_hess",
\end{verbatim}
If the option \textit{HESS}=1 is chosen in the input file, a diagonal Hessian approximation is applied for preconditioning of the gradients. If the Hessian approximation is already calculated it can be read from file (\textbf{READ\_HESS}=1). Hereby the hessian approximation is read from the folder \textit{hess}. The files written at the beginning of the current frequency stage are used, labeled for example \textit{hess/toy\_hess.vp\_200.00Hz\_it11} when the new frequency stage started with iteration 11. \\
For \textbf{READ\_HESS}=0 the diagonal Hessian approximation is calculated in IFOS3D as described in section~\ref{sec:hess}.\\
Due to runtime and storage it is recommendable to use only a subset of receivers for its calculation. The option \textbf{REC\_HESS} is not fully implemented at the moment, and it is therefore necessary to change the receiver file or the distance in the receiver array to its double or its triple value. \\
The water level (\textbf{WATER\_HESS}) used for the calculation of the preconditioning matrix (equation~\ref{Hess_precon}) can be estimated from the gradient of the first iteration and the diagonal Hessian matrix. It is not calculated directly in IFOS3D. This results in the fact, that it is not possible to perform a full inversion at a stretch, but the inversion needs to be split up into the different inversion stages.\\
The diagonal Hessian matrix is stored in the folder \textit{hess} and is labeled \\e.g. \textit{hess/toy\_hess.vp\_200.00Hz\_it11}. It is calculated at the beginning of each frequency stage and is applied within this frequency stage.
\subsection{L-BFGS}
\textcolor {Gray}{Apply\_L\_BFGS\_(LBFGS) = 0\\
Number\_of\_inverted\_parameters\_(NUMPAR) = 2\\
Number\_iterations\_used\_for\_LBFGS\_(BFGSNUM) = 5} \vspace{0.1cm}\\
\begin{verbatim}
"L-BFGS" : "comment",
"LBFGS" : "0",
"NUMPAR" : "2",
"BFGSNUM" : "5"
\end{verbatim}
If \textbf{LBFGS}=1 the L-BFGS method is applied as described in section~\ref{sec:lbfgs}. The L-BFGS approach can only be used combined with the Hessian preconditioning (\textbf{HESS}=1). So far we tested the L-BFGS approach only for the seismic velocities, not for density. The number of inverted parameter classes (\textbf{NUMPAR}) can thus be chosen as 1 ($v_p$ only) or 2 ($v_p$ and $v_s$).\\
The L-BFGS algorithm estimates the Hessian from changes in gradients and models of the previous \textbf{BFGSNUM} iterations.
......
......@@ -30,7 +30,7 @@ The model does not contain a free surface (\textbf{FREE\_SURF}=0) and as boundar
\end{center}
\end{figure}
Sources and receivers are arranged within $x$-$y$-planes, as indicated in Figure~\ref{fig:toy_model}. The sources are listed in \textit{sources/sources\_toy.dat} (\textbf{SOURCE\_FILE}). We use 12 (3$\times$4) sources in 92\,m depth with a central frequency of \textbf{FC}=300\,Hz. The sources are vertical directed point forces (\textbf{QUELLTYP}=4) with sin$^3$-wavelets (\textbf{QUELLART=4}) as source time functions. The corresponding source wavelet can be seen in figure~\ref{fig:toy_wavelet}a. This source wavelet comprises the low frequencies with $40\%$ of the maximum amplitude left at about 400\,Hz (figure~\ref{fig:toy_wavelet}b). \\
The receivers are arranged using a horizontal receiver array (\textbf{REC\_ARRAY}=1) in 24\,m depth. The distance between the receivers is \textbf{DRX}=\textbf{DRY}=10 gridpoints, which gives a total number of 169 receivers.
The receivers are arranged using a horizontal receiver array (\textbf{READREC}=2, \textbf{REC\_ARRAY}=1) in 24\,m depth. The distance between the receivers is \textbf{DRX}=\textbf{DRY}=10 gridpoints, which gives a total number of 169 receivers.
\subsection{Inversion parameters}
For the inversion we use homogeneous starting models with $v_s=3600$\,m/s and $v_p=6200$\,m/s. In total 60 iterations are performed (\textbf{ITMIN, ITMAX }= 1, 60) and five frequencies are employed simultanously (\textbf{NFMAX}=5). For preconditioning of gradients a local taper is applied around source and receiver positions (\textbf{DAMPTYPE}=2). We do not use the diagonal Hessian approximation (\textbf{HESS}=0) or the L-BFGS method (\textbf{LBFGS}=0). Out of the total number of 12 shots we use 4 shots for the steplength estimation (\textbf{NSHOTS\_STEP}=4) and start with an initial test steplength of 0.02 (\textbf{TESTSTEP}=0.02). The filenames for gradient and model output are given in the input file.
\subsubsection*{The workflow}
......@@ -119,8 +119,8 @@ Model parameters are written into the folder \textit{model} in each iteration fo
The results for the toy example is plotted for two 2D slices. Figure~\ref{fig:toy_result1} shows a horizontal slice through the box area. The three different box areas seen in the real model in a) and c) are well resolved for $v_p$ (b) and $v_s$ (d) after 60 iterations. However the $v_s$-model offers a clearly higher resolution due to the smaller wavelengths compared to the $v_p$-model. In transmission geometry a resolution down to one wavelength is possible and higher frequencies would be necessary to gain a better resolved box in $v_p$.
\begin{figure}[h!]
\begin{center}
\includegraphics[width=\textwidth]{fig_toy/toy_model_result}
\caption[Toy example - final inverted models, horizontal slice]{Final inverted models (60 iteartions) compared to real models for horizontal slice at $z$=60\,m: a) real model $v_p$, b) inverted model $v_p$, c) real model $v_s$ and d) inverted model $v_s$; }\label{fig:toy_result1}
\includegraphics[width=\textwidth]{fig_toy/toy_model_result_new}
\caption[Toy example - final inverted models, horizontal slice]{Final inverted models (60 iteartions) compared to real models for horizontal slice at $z$=60\,m: a) real model $v_p$, b) inverted model $v_p$, c) real model $v_s$ and d) inverted model $v_s$. The dashed line indicates the absorbing frame; }\label{fig:toy_result1}
\end{center}
\end{figure}
A vertical slice of the models is plotted in figure~\ref{fig:toy_result2}. Overall, the vertical direction is more difficult to reconstruct, as can be ssen in the inverted models of $v_p$ (b) and $v_s$ (d) compared to the real models (a,b). The boundaries of the box are relatively well recovered for both parameters. The inversion of $v_p$ reconstructs the two main areas of the box, however the small high velocity area (dark red) cannot be resolved. This area is indicated in the inverted $v_s$ model, again showing the higher resolution of this parameter.
......
......@@ -48,7 +48,7 @@
"RUN_MULTIPLE_SHOTS" : "1",
"Following lines can be ignored if SRCREC!=2" : "comment",
"Plane wave (PW) excitation,if PLANE_WAVE_DEPTH>0, SRCREC is treated as 0" : "comment",
"Plane wave (PW) excitation" : "comment",
"depth_of_PW_excitation_(no<=0)_(in_meter)" : "comment",
"PLANE_WAVE_DEPTH" : "0.0",
"dip_of_PW_from_vertical_(in_degrees)" : "comment",
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment