Commit 9d9866d5 authored by thomas.forbriger's avatar thomas.forbriger
Browse files

[MERGE] (master|stf_issue12): provide soutifu toy examples on master

parents 6a5d3486 6ca22f56
......@@ -70,7 +70,7 @@ documentation.
User documentation
------------------
The theory behind the Fourier domain least squares engine is outlined by
The theory behind the Fourier domain least squares procedure is outlined by
Lisa Groos (2013, Appendix F, page 146). She also describes a way to find an
approrpiate water level by application of the L-curve criterion (Groos, 2013,
Appendix G, page 148).
......@@ -87,6 +87,9 @@ tests/onlinehelp Provides access to these texts. Just issue
to get access.
A toy example and a step-by-step introduction are provided in subdirectory
src/ts/wf/testcases in README.soutifu
A short descrpition of the library and the accompanying program soutifu is
provided on the OpenTOAST web-page:
http://www.opentoast.de/Data_analysis_code_soutifu_and_libstfinv.php
......
Schicht ueber Halbraum Referenz Modell
Layer over halfspace - reference model
number of layers: 2
earth radius (km): 0.000 reference frequency (Hz): .000
......
Schicht ueber Halbraum Referenz Modell
Layer over halfspace - slightly different from reference model
number of layers: 2
earth radius (km): 0.000 reference frequency (Hz): .000
......@@ -6,8 +6,8 @@ Schicht ueber Halbraum Referenz Modell
Zb alpha beta rho Qalpha Qbeta Rb
-----------------------------------------------------------------------------
.000 .33180 .00000 .001300 1000.000 -1.000
0.0050000 0.3500 0.18000 1.500000 5.000 5.000
halfspace: 1.0000 0.55000 2.200000 20.000 20.000
0.0052000 0.4400 0.22000 1.500000 5.000 5.000
halfspace: 1.6000 0.57000 2.200000 20.000 20.000
par units description
----------------------------------------------------------------------------
......
# this is <Makefile>
# ----------------------------------------------------------------------------
# $Id$
#
# Copyright (c) 2011 by Thomas Forbriger (BFO Schiltach)
#
......@@ -34,8 +33,15 @@
all:
TESTCASEMAKE=$(filter-out %.bak,$(wildcard Makefile*))
flist: $(TESTCASEMAKE) $(wildcard *.par) $(wildcard *.tpl)
echo $^ | tr ' ' '\n' | sort > $@
READMEFILES=$(filter-out %.bak,$(wildcard README*))
flist: $(TESTCASEMAKE) $(READMEFILES) $(wildcard *.par) $(wildcard *.tpl) \
$(wildcard *.gpt)
echo $(TESTCASEMAKE) | tr ' ' '\n' | sort > $@
echo '----' >> $@
echo $(filter-out $(TESTCASEMAKE) $(READMEFILES),$^) \
| tr ' ' '\n' | sort >> $@
echo '----' >> $@
echo $(READMEFILES) | tr ' ' '\n' | sort >> $@
.PHONY: edit
edit: flist; vim $<
......@@ -52,7 +58,7 @@ clean: ;
soutifu%: Makefile.soutifu ; $(MAKE) -f $< $@
# soutifu test cases
SOUTIFUCASES=1 1_add 2a 2b 3a 3b 3c
SOUTIFUCASES=1 1_add 2a 2b 2c 3a 3b 3c
SOUTIFURESULTS=$(addprefix soutifu,$(addsuffix _report,$(SOUTIFUCASES)))
soutifutests:
$(MAKE) clean
......
# this is <Makefile.deconv>
# ----------------------------------------------------------------------------
# $Id$
#
# Copyright (c) 2011 by Thomas Forbriger (BFO Schiltach)
#
......
# this is <Makefile.filters>
# ----------------------------------------------------------------------------
# $Id$
#
# Copyright (c) 2011 by Thomas Forbriger (BFO Schiltach)
#
......
# this is <Makefile.foutra>
# ----------------------------------------------------------------------------
# $Id$
#
# Copyright (c) 2011 by Thomas Forbriger (BFO Schiltach)
#
......
# this is <Makefile.phasedsignals>
# ----------------------------------------------------------------------------
# $Id$
#
# Copyright (c) 2011 by Thomas Forbriger (BFO Schiltach)
#
......
# this is <Makefile.resaseda>
# ----------------------------------------------------------------------------
# $Id$
#
# Copyright (c) 2011 by Thomas Forbriger (BFO Schiltach)
#
......
# this is <Makefile.soutifu>
# ----------------------------------------------------------------------------
# $Id$
#
# Copyright (c) 2011 by Thomas Forbriger (BFO Schiltach)
# Copyright (c) 2011, 2016 by Thomas Forbriger (BFO Schiltach)
#
# testcases for soutifu
# Testcases for soutifu
#
# See also README.soutifu
#
# ----
# This program is free software; you can redistribute it and/or modify
......@@ -22,14 +23,19 @@
# Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
# ----
#
# Rules defined in this file require a complete installation of Seitosh along
# with Linux tools like ps2pdf, gnuplot, gs, gv, evince, sed, and a2ps
#
# REVISIONS and CHANGES
# 24/06/2011 V1.0 Thomas Forbriger
# 30/09/2011 V1.1 test convolution of additional traces
# 19/03/2016 V1.2 - provide detailed documentation
# - add noise to synthetic test data
# - provide control plots
#
# ============================================================================
#
SOUTIFUCASES=1 1_add 2a 2b 3a 3b 3c
SOUTIFUCASES=1 1_add 2a 2b 2c 3a 3b 3c
SOUTIFURESULTS=$(addprefix soutifu,$(addsuffix _report,$(SOUTIFUCASES)))
.PHONY: all
......@@ -41,38 +47,69 @@ all:
clean:
-/bin/rm -vf flist *.xxx *.xxx.* *.sff *.TZ *.TR *.grn *.su
.PHONY: soutifuclean
soutifuclean:
-/bin/rm -vf soutifu* *.sff *.su stf.fil.xxx *description.xxx
#
# produce synthetic shot gathers, using refmet
# ======================================================================
# 1. Produce synthetic shot gathers by using refmet
# =================================================
# This section provides synthetic seismograms. Model files defined in variable
# MODELS are supported. Waveforms are simulated by the reflectivity method as
# implemented in refmet.
#
MODELS=2lay.par
# 1a. Prepare control files
# -------------------------
# refmet control files are created from templates such that some parameters
# can be controlled from the command line through make variables.
#
# Basic model files are:
MODELS=2lay.par 2layb.par
# source configuration
# --------------------
# duration of source wavelet
SRCDUR=0.01
# depth of source
SRCDEPTH=0.0005
# control file for single vertical force
zforce.xxx: zforce.tpl
sed -e 's/<SDUR>/$(SRCDUR)/' \
-e 's/<SDEPTH>/$(SRCDEPTH)/' $< > $@
# main file
# ---------
# Create main control file
# ------------------------
# Three main control files are provided
REFMETSRC=zforce.xxx
REFMETMOD=2lay.par
REFMETRCV=linrcv.par
.PRECIOUS: $(REFMETSRC)
#
# a) Allows all templates to be replaced
# This is not used for the testcases below
refmetmain.xxx: refmetmain.tpl $(REFMETMOD) $(REFMETSRC) $(REFMETRCV)
sed -e 's/<MODEL>/$(REFMETMOD)/' \
-e 's/<SOURCE>/$(REFMETSRC)/' \
-e 's/<RECEIVERS>/$(REFMETRCV)/' $< > $@
#
# b) Control file with standard vertical force an reference model
2laymain.xxx: refmetmain.tpl 2lay.par $(REFMETSRC) $(REFMETRCV)
sed -e 's/<MODEL>/2lay.par/' \
-e 's/<SOURCE>/$(REFMETSRC)/' \
-e 's/<RECEIVERS>/$(REFMETRCV)/' $< > $@
#
# c) Control file with standard vertical force an modified reference model
2laybmain.xxx: refmetmain.tpl 2layb.par $(REFMETSRC) $(REFMETRCV)
sed -e 's/<MODEL>/2layb.par/' \
-e 's/<SOURCE>/$(REFMETSRC)/' \
-e 's/<RECEIVERS>/$(REFMETRCV)/' $< > $@
# Simulate seismograms
# --------------------
# Based on the control files defined above, waveform data is computed by
# refmet. The files 2lay.TZ and 2layb.TZ are used as a basis for testcases
# defined below.
.PRECIOUS: seismo.TZ seismo.TR
seismo.TZ seismo.TR: refmetmain.xxx
refmet -v -4 -o seismo -c $<
......@@ -83,10 +120,66 @@ seismo.TZ seismo.TR: refmetmain.xxx
2layb.TZ 2layb.TR: 2laybmain.xxx
refmet -v -4 -o 2layb -c $<
#----------------------------------------------------------------------
# Plot waveform data
# ------------------
# Create a graphical display of waveform data used for testcases.
soutifubasewaveform.ps: 2lay.TZ 2layb.TZ
refract -d $@/cps -C -Tt "Basic synthetic time series" \
-Sm 2 -Se 2.3 -Sc 5. -Sa 3.5 -Eg -St 0.,1. $^
# ======================================================================
# 2. Fourier Bessel control operation
# ===================================
# This section provide Fourier-Bessel expansion coefficients for the suburface
# models and the vertical single force defined above. This is provided as an
# alternative display of wave propagation properties as resulting from the
# subsurface structure defined by refmet models.
#
# Use the following commands to create Postscript files:
#
# Fourier-Bessel expansion coefficients for models 2lay and 2layb:
#
# make -f Makefile.soutifu soutifu2laysyggrn.ps soutifu2laybsyggrn.ps
#
# Difference in Fourier-Bessel expansion coeffcients for both models:
#
# make -f Makefile.soutifu soutifudiffgrn.ps
# produce seismogram file by windowing such that backpropagating
# cut-off phase is removed
# syg
# ---
# Calculate Fourier-Bessel expansion coefficients
$(patsubst %.par,%syg.grn,$(MODELS)): %syg.grn: %.par zforce.xxx
syg $^ $@ -P 100,100,100.,10.
# Calculate difference between coefficients for both basic models
diff.grn: 2laysyg.grn 2laybsyg.grn; rm -fv $@; gredi $^ $@
# grepg
# -----
# Create Postscript plot of Fourier-Bessel coefficients
soutifu%grn.ps: %.grn
grepg $< -W -d $@/ps -p 0.8 -s -f 10 -i -g
# ======================================================================
# 3. Preprocess waveform data
# ===========================
#
# This section provides preprocessed waveform data as a basis for tests to be
# applied to soutifu.
#
# 3a. Remove numerical cut-off phase
# ----------------------------------
# The reflectivity code produces a backpropagating cut-off as a numerical
# artefact. This is removed by selecting a given number of samples at the
# beginning of the seismograms. As a result the following waveform data files
# are provided:
#
# synthetic2lay.sff
# synthetic2layb.sff
# syntheticseismo.sff
#
# The number of samples to be used is defined by make variable WINLIM
#
WINLIM=2000
win.fil.xxx: ; printf "win 1,$(WINLIM)\nend\n" > $@
newsynthetic%.sff: %.TZ win.fil.xxx
......@@ -97,59 +190,116 @@ synthetic2lay.sff synthetic2layb.sff syntheticseismo.sff:
$(MAKE) new$@ -f Makefile.soutifu
/bin/cp -vpd new$@ $@
# apply a filter, simulating a source time function
# 3b. Apply a filter, simulating a source time function
# -----------------------------------------------------
# As a second step to simulate recorded data, a filter is applied to the
# waveform data, which simulates a source time function or the response of a
# recording system or both. The filter control file is later used as a
# reference to study the capability to reproduce the filter characteristic by
# soutifu. The filter characteristic can be controlled by make variables.
#
# Low-pass period
STFPER=0.04
STFFAC=0.1
stf.fil.xxx: ; printf "lpb $(STFPER),4\ndel 0.03\nfac $(STFFAC)\nend\n" > $@
data%.sff: synthetic%.sff stf.fil.xxx
# Scaling factor
STFFAC=2.5
# Trigger advance
STFDELAY=0.03
stf.fil.xxx: ; printf "lpb $(STFPER),4\ndel $(STFDELAY)\nfac $(STFFAC)\nend\n" > $@
filtered%.sff: synthetic%.sff stf.fil.xxx
tidofi -v -o -cf $(word 2,$^) $@ $<
#----------------------------------------------------------------------
# additional series to be convolved
# 3c. Add noise
# -------------
# As a first step to simulate recorded data, white Gaussian noise can be added
# to the synthetic waveforms.
#
# The rms-amplitude of noise is specified by variable NOISEAMP
# a reasonable value to start with is NOISEAMP=1.e-5
NOISEAMP=0.
data%.sff: filtered%.sff
printf "noi $(NOISEAMP)\nend\n" | tidofi -v -o -cs $@ $<
# Display waveforms
# -----------------
# Display waveform data as provided for tests. This way different settings for
# the source time function simulation and data noise can be tested.
#
# soutifusyntheticplot.ps
# soutifunoisyplot.ps
# soutifudataplot.ps
soutifu%plot.ps: %2lay.sff %2layb.sff
refract -d $@/cps -C -Tt "Synthetic time series" \
-Sm 2 -Se 2.3 -Sc 5. -Sa 3.5 -Eg -St 0.,1. $^
# Display waveform amplitudes
# ---------------------------
# Output rms and peak-to-peak amplitudes. This is used to provide a reference
# to adjust the noise amplitude parameter.
%.amp.xxx: %.sff
sigval -format "%F (trace #%NT at %OFF): rms= %RMS P-P= %PPA" \
$< > $@
%.ampdat.xxx: %.sff
sigval -format "%NT %RMS %PPA" $< > $@
%.amp: %.amp.xxx; cat $<
soutifuamp.ps: plotsoutifuamp.gpt \
synthetic2lay.ampdat.xxx synthetic2layb.ampdat.xxx \
data2lay.ampdat.xxx data2layb.ampdat.xxx
gnuplot $<
# ======================================================================
# 4. Additional series to be convolved
# ====================================
# Create additional test time series which will be used to test soutifus
# ability to reproduce the actual source time function.
additional.sff: win.fil.xxx
/bin/rm -fv $@
tesiff -r 1000. -n $(WINLIM) $@
# Apply the source time function filter as defined above to the additional
# test waveform data. The filter result will be compared fo time series
# filtered by soutifu.
additionalfiltered.sff: additional.sff stf.fil.xxx
tidofi -v -o -cf $(word 2,$^) $@ $<
#----------------------------------------------------------------------
# operate on SeismicUn*x data
# ======================================================================
# 5. Operate on SeismicUn*x data
# ==============================
# In order to provide test data and test results in Seismic Un*x, such that
# they can be further processed by SU tools, we provide conversion rules.
%.su: %.sff; any2any --verbose --itype sff --otype su --overwrite $@ $<
%_conv.sff: %.su; any2any --verbose --otype sff --itype su --overwrite $@ $<
#----------------------------------------------------------------------
# Fourier Bessel control operation
# syg
# ---
$(patsubst %.par,%syg.grn,$(MODELS)): %syg.grn: %.par zforce.xxx
syg $^ $@ -P 100,100,100.,10.
# grepg
# -----
%grn.ps: %.grn
grepg $< -W -d $@/ps -p 0.8 -s -f 10 -i
# ============================================================================
# test cases
# ----------
# 6. Test cases
# =============
#
# 6a. Test case 1
# ---------------
define TESTCASE1
Synthetics are calculated with a standard source time function which is band
limited and has a duration of $(SRCDUR) seconds. The recorded data is
simulated from this by applying a low-pass filter, a scaling factor, and a
time delay. Noise with an rms-amplitude of $(NOISEAMP) additionally is added
to the simulated recordings. Because simulated recordings and synthetics have
the same wave propagation characteristics, the simulated field data can
exactly be reproduced by the source time function correction filter, at least
in the case where no noise is added to the data.
endef
export TESTCASE1
# test case 1
# -----------
SOUTIFU1SHIFT=0.4
SOUTIFU1=fdlsq:tshift=$(SOUTIFU1SHIFT)
soutifu1_data.sff: data2lay.sff; /bin/cp -vpd $< $@
soutifu1_synthetic.sff: synthetic2lay.sff; /bin/cp -vpd $< $@
soutifu1_signaldesc.xxx: stf.fil.xxx
echo "synthetics are calculated with a standard source" > $@
echo "time function which is band limited and has a" >> $@
echo "duration of $(SRCDUR) seconds. The recorded data" >> $@
echo "is simulated from this by applying low-pass filter" >> $@
echo "and a time delay:" >> $@
soutifu1_signaldesc.xxx: stf.fil.xxx zforce.xxx
printf "$$TESTCASE1" > $@
printf "\n\nContents of $< are:\n" >> $@
cat $< >> $@
printf "\n\nContents of $(word 2,$^) are:\n" >> $@
cat $(word 2,$^) >> $@
soutifu1_convseis.su soutifu1_stf.su soutifu1_add.su: \
soutifu1_data.su soutifu1_synthetic.su additional.su
......@@ -194,7 +344,7 @@ soutifu1_add_%.ps: soutifu1_add.su additional.su additionalfiltered.su
soutifu1_stf.ps: soutifu1_stf.su
stuplox -ty su -d $@/ps -E -m 0. -Y 0.9 -S -$(SOUTIFU1SHIFT) \
-R 0.6 -L 0.1 -T "soutifu test case 1: $(SOUTIFU1)" -l 3,3,3\
-R 0.7 -T "soutifu test case 1: $(SOUTIFU1)" -l 3,3,3\
-u "source time function" $<
#----------------------------------------------------------------------
......@@ -210,8 +360,21 @@ soutifu1_add_report.pdf: soutifu1_description.ps soutifu1_section.ps\
gs -sDEVICE=pdfwrite -dNOPAUSE -dBATCH -sOutputFile=$@ $^
#----------------------------------------------------------------------
# test case 2
# -----------
# 6b. Test case 2
# ---------------
define TESTCASE2
We take the same data like in test case 1. This time, however, data and
synthetics are swapped. Since synthetics are narrow band with respect to
data an appropriate water level is required.
As a consequence of swapping the role of data and synthetics, the source time
function correction necessarily is entirely acausal. Since simulated data and
synthetics are both based on the same wave propagation characteristics,
soutifu should be able to almost reproduce the simulated data, at least in the
case where no noise is added.
endef
export TESTCASE2
SOUTIFU2SHIFT=0.4
SOUTIFU2WL=1.
SOUTIFU2=fdlsq:tshift=$(SOUTIFU2SHIFT):waterlevel=$(SOUTIFU2WL)
......@@ -224,14 +387,9 @@ soutifu2x_convseis.su soutifu2x_stf.su: \
--type su $(SOUTIFU2) $^
soutifu2x_description.xxx: soutifu1_signaldesc.xxx
echo "soutifu test case 2: $(SOUTIFU2)" > $@
echo >> $@
echo "We take the same data like in test case 1." >> $@
echo "This time, however, data and synthetics are exchanged." >> $@
echo "Since synthetics are narrow band with respect to data" >> $@
echo "an appropriate water level is required." >> $@
echo >> $@
echo "Signal description from test case 1:" >> $@
printf "soutifu test case 2: $(SOUTIFU2)\n\n" > $@
printf "$$TESTCASE2" >> $@
printf "\n\nSignal description from test case 1:\n" >> $@
cat $< >> $@
soutifu2x_section.ps: soutifu2x_data.su soutifu2x_synthetic.su \
......@@ -262,8 +420,25 @@ soutifu2c_report.pdf:
/bin/mv -fv soutifu2x_report.pdf $@
#----------------------------------------------------------------------
# test case 3
# -----------
# 6c. Test case 3
# ---------------
define TESTCASE3
The model for which synthetics are calculated produced a slightly smaller wave
propagation. It is no longer possible to repoduce the data with a single
source time function correction filter for all traces. Therefore a compromise
has to be made between far and near offset traces. This test case can be used
to study the effect of waterlevel and offset dependent weights.
Signal description:
Synthetics are calculated with a standard source time function which is band
limited and has a duration of $(SRCDUR) seconds. The recorded data is
simulated from this by applying a low-pass filter, a scaling factor, and a
time delay.
Noise with an rms-amplitude of $(NOISEAMP) is added to the data.
endef
export TESTCASE3
SOUTIFU3SHIFT=0.4
SOUTIFU3WL=1.e-4
SOUTIFU3EXP=0.
......@@ -276,22 +451,13 @@ soutifu3x_convseis.su soutifu3x_stf.su: \
soutifu -v -o -wc soutifu3x_convseis.su -ws soutifu3x_stf.su \
--type su $(SOUTIFU3) $^
soutifu3x_description.xxx: stf.fil.xxx
echo "soutifu test case 3: $(SOUTIFU3)" > $@
echo >> $@
echo "The model for which synthetics are calculated is too" >> $@
echo "slow in its wave propagation properties. Therefore" >> $@
echo "a compromise has to be made between far and near" >> $@
echo "offset traces." >> $@
echo >> $@
echo "Signal description:" >> $@
echo "Synthetics are calculated with a standard source" >> $@
echo "time function which is band limited and has a" >> $@
echo "duration of $(SRCDUR) seconds. The recorded data" >> $@
echo "is simulated from this by applying low-pass filter" >> $@
echo "and a time delay. The corner period of the low-pass is" >> $@
echo "so small that it pracically has no effect:" >> $@
soutifu3x_description.xxx: stf.fil.xxx zforce.xxx
printf "soutifu test case 3: $(SOUTIFU3)\n\n" > $@
printf "$$TESTCASE3" >> $@
printf "\n\nContents of $< are:\n" >> $@
cat $< >> $@
printf "\n\nContents of $(word 2,$^) are:\n" >> $@
cat $(word 2,$^) >> $@
soutifu3x_section.ps: soutifu3x_data.su soutifu3x_synthetic.su \
soutifu3x_convseis.su
......@@ -321,15 +487,24 @@ soutifu3c_report.pdf:
/bin/mv -fv soutifu3x_report.pdf $@
#======================================================================
# 7. Compile reports
# ==================
#
%_description.ps: %_description.xxx
a2ps -o $@ -r --center-title=$(patsubst %_description.xxx,%,$<) \
-l 60 $<
a2ps -o $@ --center-title=$(patsubst %_description.xxx,%,$<) $<
%_report.ps: %_description.ps %_section.ps %_stf.ps
gs -sDEVICE=pswrite -dNOPAUSE -dBATCH -sOutputFile=$@ $^
%_report.pdf: %_description.ps %_section.ps %_stf.ps
gs -sDEVICE=pdfwrite -dNOPAUSE -dBATCH -sOutputFile=$@ $^
# ======================================================================
# 8. Previewer
# ============
#
PDFPREVIEWER=evince
%.pdf: %.ps; ps2pdf $<
%.psp: %.ps; gv $<; /bin/rm -fv $<
%.pdp: %.pdf; acroread $<; /bin/rm -fv $<
%.pdp: %.pdf; $(PDFPREVIEWER) $<; /bin/rm -fv $<
# ----- END OF Makefile.soutifu -----
# this is <Makefile.teseco>
# ----------------------------------------------------------------------------
# $Id$
#
# Copyright (c) 2011 by Thomas Forbriger (BFO Schiltach)
#
......
this is <README.soutifu>
============================================================================
soutifu toy example - step by step
==================================
Purpose and basic usage
-----------------------
src/ts/wf/testcases/Makefile.soutifu contains the definition of test cases for
soutifu. Change directory to src/ts/wf/testcases and execute
make soutifutests
This will run all tests as defined in Makefile.soutifu The results will be
presented in PDF files:
soutifu1_add_report.pdf soutifu2b_report.pdf soutifu3b_report.pdf
soutifu1_report.pdf soutifu2c_report.pdf soutifu3c_report.pdf
soutifu2a_report.pdf soutifu3a_report.pdf
These reports contain descriptive texts and waveform displays. They serve as a
quick test for proper operation and should be useful as a first approach to
soutifu operation. The setup of the test cases can be controlled through make
variables passed on the command line. This way they can serve as toy examples.
Principle of operation for test cases and toy examples
------------------------------------------------------
The test cases are applied to simulated recordings and corresponding
synthetics. Both are computed with the reflectivity algorithm as implemented
in the program refmet.
A plot of the waveform data all other components are based on, can be produced
by
make soutifubasewaveform.pdf
To study the surface wave propagation properties, plots of Fourier-Bessel
expansion coefficients are provided too. Additionally a plot of the difference
between expansion coefficients for both subsurface models (2lay and 2layb) ㄞs
provided:
make soutifu2laysyggrn.pdf soutifu2laybsyggrn.pdf
make soutifudiffgrn.pdf