HOW.TO-gcf4ag

# GLOSSARY
# sv - seismic velocity
# sax - seismic acceleration
#

PURPOSE: Determine transfer functions and Wiener filters between a (3-D) seismometer series and an AG drop residual

# You need a combined seismogram:cd Seismo/gcfsource append-seismograms-mc.env# rm d/tmp.mc# foreach c ( Z E N )# tslist-app -O:`label SV,${c}` d/tmp.mc -I + \# G2011.166/GCF.2011:166:00:[0-5]0:00.3U93${c}2.ts \# G2011.166/GCF.2011:166:0[1-9]:*.3U93${c}2.ts \# G2011.166/GCF.2011:166:1[0-1]:*.3U93${c}2.ts \# G2011.166/GCF.2011:166:12:00:00.3U93${c}2.ts# end
# Actually, we should phase-correct the seismograms
# This is done in the version of Jan. 2013 of append-seismograms-mc.env

# You need a drop residual. Cut out a suitable piece
tslist AG/o/scg-cal-drops.ra.ts -qqq \
       -BHc2011,6,15,0,20,00 -C4 -Ff10.3 -U2011,6,15,12,09,05 \
       -w AG/scg-cal-drops.ra.tsf

# You need a tsfedit-cont file: AG/cmd.all.drop.times.tsi
#
# ${TSECOMMAND} At 2011 06 11 13 09 34 000
#
# from which a piece suitable for the period will we extracted

# First, prepare a -1024 ... +1024 suite of shifted seismograms
# resampled at the AG drop times: First vertical accelerations, then
# the two horizontal components. (We might add a horizontal-speed file)

sample-gcf4AG -D sax -E sample-gcf4AG.tse,DIFF \
              'G2011.166/GCF.2011:166:??:??:00.3U93Z2.ts'

sample-gcf4AG -D sv 'G2011.166/GCF.2011:166:??:??:00.3U93N2.ts'
sample-gcf4AG -D sv 'G2011.166/GCF.2011:166:??:??:00.3U93E2.ts'

# HOW To use sample-gcf4AG? Do   m sample-gcf4AG
# Here is an example for producing 2Hz sax data, 50 lags in every direction
sample-gcf4AG -D 2Hzsax -I 50 -l -50,51 \
      -E sample-gcf4AG.tse,DIFF \
      -E sample-gcf4AG.tse,LP2HZ \
      'G2011.166/GCF.2011:166:??:??:00.3U93Z2.ts'

# Then we compute covariance, power spectrum, and cross-spectrum,
# the latter will create a Wiener filter bank.
#
# The seismogram's power spectrum:
# A difference filter -.99 might be needed, perhaps even greater

cd ~/Seismo/gcf
sasm03 \!sasm03.ins

# (edit sasm03.ins to activate the DYDT filter for sax)
# In the first show of the Bartlett Power spectrum, save the
# spectrum to file (e.g. the sv east component with W svE.psp )
# Copy the power spectrum (should have domain width 1024) to AG/

echo " dB " >! AG/svE-1024.psp
cat svE.psp >> AG/svE-1024.psp

# Covariance

# The programs that calculate covariance and correlation are
#    sax-ag-corr     sax-sax-corr     sv-ag-corr     sv-sax-corr     sv-sv-corr
#
# In sv-ag-corr change set comps = ( x ) appropriately, e.g. x = E
# In sasm06-sv-ag.ins   enter component E mutatis mutandis
#
# The task consists of computing covariance series and, together with the power spectrum computed before (by sasm03),
# compute transfer spectra and Wiener filters. We hope that the correlations increase when we run the seismometer
# records through the Wiener filters.

cd ~/Seismo/gcf/AG
set corrdata=1
source sv-ag-corr

# For the Z-component we use acceleration, sax-ag-corr and sasm06-sax-ag.ins, mutatis mutandis

sasm06 \!sasm06-sv-ag.ins

# where at the stage of Impulse response the Wiener filter bank
# ( svE-ag.wfb ) will be output. Apply the Wiener filters:
#

sample-gcf4AG -T agp -P sax -O AG/saxdata -l -161,80 \
              -E sample-gcf4AG.tse,DIFF                \
              -E sample-gcf4AG.tse,WFSAX                \
              'G2011.166/GCF.2011:166:??:??:00.3U93Z2.ts'

sample-gcf4AG -T agp -P sv -O AG/svdata -l -161,80    \
              -E sample-gcf4AG.tse,WFSV                 \
              'G2011.166/GCF.2011:166:??:??:00.3U93E2.ts'

sample-gcf4AG -a -T agp -P sv -O AG/svdata -l -161,80 \
              -E sample-gcf4AG.tse,WFSV                 \
              'G2011.166/GCF.2011:166:??:??:00.3U93N2.ts'
#
#OBS!!! ________________________________
# Third call runs option -a in order to append the   |
# mc-files that are created.                                     |
#---------------------------------------------------
#
# For the 2Hz variant you will need

sample-gcf4AG -T agp -P 2Hzsax -I 50 -O AG/2Hzsaxdata -l -25,35 \
              -E sample-gcf4AG.tse,DIFF         \
          -E sample-gcf4AG.tse,LP2Hz         \
              -E sample-gcf4AG.tse,WFSAX          \
              'G2011.166/GCF.2011:166:??:??:00.3U93Z2.ts'

sample-gcf4AG -T agp -P 2Hzsv -I 50 -O AG/2Hzsvdata -l -25,25 \
          -E sample-gcf4AG.tse,LP2Hz         \
              -E sample-gcf4AG.tse,WFSV           \
              'G2011.166/GCF.2011:166:??:??:00.3U93E2.ts'

sample-gcf4AG -a -T agp -P 2Hzsv -I 50 -O AG/2Hzsvdata -l -25,25 \
          -E sample-gcf4AG.tse,LP2Hz         \
              -E sample-gcf4AG.tse,WFSV           \
              'G2011.166/GCF.2011:166:??:??:00.3U93N2.ts'

Drop scatter reduction

# Now you can produce a plot showing the drop scatter (rms) reduction
# for the different components

source ag-sv-wf-predregres

# Edit the component parameter,
# also edit the plot ranges, ticks etc. in this script!
#
# Plots are found like this one: plot/ag-saxZ-wf-predres-coeff.ps
# It shows regression the coefficients (amplitude adjustment of the
# Wiener-predicted time series) and the RMS difference after
# this simple regression exercise.

# Find out the best correlated seismograph series (what particular shift?)
# Answer:
awk '{print $1,$2,$11-$20}' ag-saxZ-wf-predres.rms | sort -nr -k3 | m
awk '{print $1,$2,$11-$20}' ag-svE-wf-predres.rms | sort -nr -k3 | m
awk '{print $1,$2,$11-$20}' ag-svN-wf-predres.rms | sort -nr -k3 | m

# Take the series from the top of these lists into urtap-wfp.ins
# and compute regression. Iterate for outliers if you wish

urtapt \!urtap-wfp.ins
#
# RESULTS: ADMITTANCES, LOCAL COPHASES
#
# SITE/FILE: agp-s-004.mc lon/lat:    11.9260   57.3964
#
# Normalized Chi^2 of fit : 1.88D+02, R=   1.37D+01 , X_1 X_2 = 1.00 2.30, Nev= 5
#
# The error information below is compatible with a unit normalized Chi^2.
# Gain factor appears rectified for the NDR
# <rslist>d> rmsu= 13.6934412298610 weff= 1.00000000000000 rchi= 1.00064262960750
#
#      #b Dominating tide   Frequ.    Amplit.   Phase     Co..admittance parameter..Quad +- 68.3% conf   Gain   Cophase
#          argum.numbers     [cyc/d]    [m]      [deg]
# SVN   1 <svdata/agp-s+013.mc>                             -0.182995 +-0.037891
# SVE   2 <svdata/agp-s+000.mc>                             -0.272778 +-0.036064
# SAXZ 3 <saxdata/agp-s-004.mc>                            -0.447830 +-0.035368
# /     4 Linear [m]/(8510*1.38889D-03[h])                  -6.479847 +-1.177031
# /     4 Linear [m]/[year] ~/[h]                         -4.8058D+03 +-8.73D+02 -5.4824D-01 +-9.96D-02
# --    5 Const 1.000D+00                                   0.003429 +-0.280164
#

# What's the drop scatter reduction (RMS)?
# Before fit:

tslist saxdata/agp-s+000.mc -L'A|R' -I -Erms.tse,R | fgrep RMS-dev

# After fit:
tslist o/ag-wfp.ra.ts -I -Erms.tse,R | fgrep RMS-dev

# Not too bad, ey?

CYCLE SLIPS

# There's the suspicion that AG mesurements are subject to cycle slips
# see and source ~/Seismo/AG/optimum_cycslip
# Find a cycle height from a compromise of good RMS-reduction and
# correlation-increase in the plot plot/cycls*.ps
# cf HOW.TO-cycleslip

# (#)
set cyci=8.
set lag=-004
setenv COMP Z
tslist ${SAXD}data/${SAXD}-s${lag}.mc -L'A|R' -L'S|'$COMP -I \
       -CYCSLIP${cyci},-2,2,O#51:cycsl-${COMP}.dat \
       -O:`label AG,CYC${COMP}` ${SAXD}data/ag-dropres-cyc.mc | m

# This outputs the cycleslips in ASCII to ./cycsl.dat
# and the cycslip-adjusted AG data to ./saxdata/ag-dropres-cyc.mc
# (.mc because we will produce a range of different cycslip solutions).

# This solution (Z H) uses SAX and horizontal H:

tslist tmp.mc -L'A|R' -L'S|Z' -L'S|H' -I \
       -CYCSLIP11.0,-2,2,T=0. -O:`label AG,CYCZRR` saxdata/ag-dropres-cyc.mc

# ---
# For the Z H prediction:
# add option      -HR          in sample-gcf4ag
# change label to 'AG|CYCZRR'       "      "
# OBS! The label for horizontal speed in sv-ag-corr must be 'R'
# ---

tslist tmp.mc -L'A|R' -L'S|Z' -L'S|E' -L'S|N' -I \
       -CYCSLIP10.0,-2,2,+S,T=0.,O#51:cycsl-ZEN.dat \
       -O:`label AG,CYCZEN` saxdata/ag-dropres-cyc.mc
#
# edit sax-ag-corr
set getcorrdata=1
source sax-ag-corr
# and assign setenv AGRES CYCZEN
source sv-ag-corr
sasm06 \!sasm06-sax-ag.ins

# 2011-AUG-24: The above shows too high coherence above 45 Hz, since the sax-power spectrum has a dip there.
# Actually, make sure you process sax and activate TSF EDIT DYDT! The dip was for N-velocity
# 47 .. 49 Hz: 40 dB,
# 49 .. 50 Hz: flat
# Need to prepare a better one with sasm03
# It's also worth a try at this point to do multi-channel linear prediction ( ~/sas/p/mt/linpred.f )
#
# # Interlude, again sasm03
# # In the first show of the Bartlett Power spectrum, save the
# # spectrum to file (e.g. the sv east component with W sax.psp )
# # Copy the power spectrum (should have domain width 1024) to AG/
#
# echo " dB " >! AG/sax-1024.psp
# cat sax.psp >> AG/sax-1024.psp
#
# 2011-AUG-24: This time coherence became really sparse. We've obtained a Wiener filter bank
# sax-ag.wfb

cd ../
sample-gcf4AG -AG AG/saxdata/ag-dropres-cyc.mc 'AG|CYCZEN' -T agp -P sax \
              -O AG/saxdata -l -161,80 \
              -E sample-gcf4AG.tse,DIFF \
              -E sample-gcf4AG.tse,WFSAX \
              'G2011.166/GCF.2011:166:??:??:00.3U93Z2.ts'
# using the -AG option to copy the cycslip-reduced drop-residual into the resulting .mc files.
# Creates AG/saxdata/agp-s$toff.mc
# 2011-AUG-24: O.K. so far.

# source ag-sv-wf-predregres   no, rather:
source ag-sax-wf-predregres

sample-gcf4AG -AG AG/saxdata/ag-dropres-cyc.mc 'AG|CYCZH' -T agp -P sv \
              -O AG/svdata -l -161,80   \
              -E sample-gcf4AG.tse,WFSV \
              'G2011.166/GCF.2011:166:??:??:00.3U93H2.ts'

# Edit urtap-wfp.ins and choose the drop residual to to fit to. E.g.
# 21 ^ saxdata/agp-s+000.mc fmt='L:A|R' (alright with any agp-s+???)
# 21 ^ o/ag-dropres.ts      fmt='U'

urtapt \!urtap-wfp.ins
source resid-histograms

# Satisfied? If in doubt, try different -CYCSLIP options in tslist at (#)