Introduction

The purpose of this tutorial is to generate surface-wave noise for a known 1-D structural model by applying randomly oriented point sources at the surface of the structure and then propagating the 3-D ground motion to two receivers. The ground motions consist of multimode Rayleigh and Love waves.

After generating a realistic noise sequence of 900 seconds length with a sampling interval of 0.002 sec, these noise segments are then cross-correlated and stacked to form the empirical Green’s function. The noise generation is done by the DOITHF script and the cross-correlation and stacking by the DOSTKZNE script.

This tutorial is set up to emulate real data acquisition. The total record of 900 sec at 0.002 sec sampling and receivers 100 meters apart  is reasonable for a field study.

Caution

Many details of processing to determine the empirical Greens functions by cross-correlation have been simplified. For real data I have only applied a frequency domain whitening before cross-correlation. Since tis simulation assumes noise sources equally spaced in origin time, the records may be very spiky. I had no luck with these simulations. However, if I applied a AGC to the cross-correlation time series before stacking, then I got something that gave acceptable results for analysis using do_mft.   The value of these scripts is that then be modified through the use of different velocity models or different distances between the two receivers.  As the seapration distance is made smaller, I would expect that it will be difficult to get dispersion data at the longer periods. There may also be a problem with the many modes arriving, so that a dispersion mode may not be able to be identified.

The other caution is that this script is set up to emulate a high frequency data acquisition where cartesian coordinates are appropriate. Anotehr script will have to be written to emulate data sets for stations widely separated, e.g., the TArray of USArray, where geocentric coodrinates have to be used.

Download and distribution

This exercise is contained in the gzip’d tar file HFAMBIENTNOISE.tgz

After downloading execute the following commands:
HFAMBIENTNOISE.dist/
HFAMBIENTNOISE.dist/DISP.PLT
HFAMBIENTNOISE.dist/DOITHF
HFAMBIENTNOISE.dist/DOSTKZNE
HFAMBIENTNOISE.dist/Models/
HFAMBIENTNOISE.dist/Models/CUS.mod
HFAMBIENTNOISE.dist/Models/soil.mod
HFAMBIENTNOISE.dist/Models/soilm2.mod
HFAMBIENTNOISE.dist/Models/tak135sph.mod


DOITHF

This shell script performs the simulation. It has many comment lines to describe each operation. The initial part of the shell script defines the parameters that control the results. The reason for writing the script in this manner is so that much of it can be used to simulate the noise fields at two stations so that ambient noise cross-correlation techniques could be used to estimate the empirical Green's function.

For the example here, we consider the velocity model HFAMBIENTNOISE.dist/Models/soilm2.mod. The dispersion curves for this model are given Figure 1. The fundamental and higher modes are plotted. In this simulation the effects of Q are not considered since the eigenfunction programs were executed tith the -NOQ flag, e.g., sregn96 -NOQ and slegn96 -NOQ.  The group velocities indicate that there could be arrivals propagating as slow as 80 m/s. (0.08 km/sec).

disp.png
Fig. 1. Dispersion curves for the test velocity model.  The first row gives phase velocity and the second gives group velocity. The left column is for Love waves and the right is for Rayleigh waves.


As set up, the DOITHF will generate 900 sec of noise for each of two stations.

The script DOITHF is built on the DOITHVHF described under other tutorials. The only difference is that the creation of the synthetic is placed into a function, and noise is simulated at two stations. The output files have names such as 1.010101.Z and 2.010101.Z  which are the vertical component noise simulations for stations 1 and 2, respectively for time segment 010101.  All segments are then stacked to form the 900 second long noise waveform 1.Z.stk and 2.Z.stk

The initial pare of DOITHF  sets up the parameters for the simulation:

#!/bin/bash

#####
# create a noise data set at two stations by generating synthetic
# motions from randomly distributed point forces
# applied at the surface
#
# To mimic actual field recordings, a long time series is
# created which will then be analyzed by a separate script.
#
# This script is designed for high frequency
# motions for local site studies
#####

#####
# Velocity model
#####
VMOD="Models/soilm2.mod"
#####
# for surface wave synthetics
#####
NMODE=100

#####
# noise sources: These occur in the region
# -XMAX <= x <= XMAX
# -YMAX <= y <= YMAX
# except about a region DMIN about the receiver
#####
DELTA=0.002 # sample interval in seconds

DMIN=0.025 # exclude sources within DMIN km of receiver
XMAX=1.000 # define source region
YMAX=1.000

#####
# define the receiver coordinate By placing the receivers along the x-axis, the Rayleigh wave fill be on the Z and E components
# XR,YR and the Love on the N component after the cross-correlations and stacking
#####
XR1=0.050
YR1=0.0
XR2=-0.050
YR2=0.0

#####
TMAX=900 # The sac file will 0 to TMAX seconds long
#
NSRC=10000 # number of random sources
# these will occur at intervals of TMAX/NSRC seconds

The next part of the script defines some functions. The details of the functions are not given in this description.

function getvmodextreme () {
#####
# examine the velocity model to determine the
# minimum and maximum shear velocities which
# will be used for the noise sampling
#
# the following are returned globally:
# VMIN - minimum S velocity in the model
# VMAX - maximum S velocity in the model
#####
}

function getsrc()
{
#####
# get random coordinates in the region
# -XMAX <= x <= XMAX
# -YMAX <= y <= YMAX
#
# the following are returned globally:
# (XS,YS) - source coordinates in km
# (EVAL,EVLO) - sourc coordinate in geocentrc coordinates
# - for simplicity the receivers are assumed to be
# - near (0,0) so that te conversion from km to degree
# - does is essentially cartesian
#####
}


function getforce()
{
#####
# get the components for the force to be applied at the surface
# the following are returned globally:
# FN - force component in north direction
# FE - force component in east direction
# FD - force component in down direction
#
#####
}

function getdistaz ()
{
#####
# get distance, azimuth, backazimuth
# Input arguments:
# 1 Source x coordinate
# 2 Source y coordinate
# 3 Receiver x coordinate
# 4 Receiver y coordinate
#
# returns the following global variables
# DIST - distance between source and receiver in km
# AZ - azimuth from source to receiver in degrees
# BAZ - backazimuth from receiver to source in degrees
#####
}

function makesyn ()
{
#####
# make a synthetic for the global location and global force
# for a particular distance, azimuth, backazimuth
#####
# input:
# $1 = distance
# $2 = azimuth
# $3 = backazimuth
# $4,$5 x,y coordinates of source
# $6,$7 x,y coordinates of receiver
#
# return
# sac files T.Z T.N T.E
#####
}

The last part of the script performs the simulation:

##### everything below here does the synthetic of the noise #####

#####
#    clean up previouse run
#####
rm -f ??????.stk
rm -f *.Z
rm -f *.N
rm -f *.E

#####
#    get the extreme values of the S velocity from the model
#####
getvmodextreme
echo VMIN=$VMIN VMAX=$VMAX

#####
#    first generate the eigenfunctions so that the
#    synthetics can be made
#    The time window must be long enough to encompass the
#    arrivals at the fastest and slowest velocities
#####
NPT=`echo $XMAX $YMAX $VMIN $VMAX $DELTA | awk \
        '{ DIST=sqrt($1*$1 + $2*$2) ; TWIN=(DIST/$3 ) ; print int(TWIN/$5)}' `
DIST=`echo $XMAX $YMAX | awk '{print sqrt($1*$1 + $2*$2)}' `
echo DIST $DIST XMAX $XMAX YMAX $YMAX
cat > ddfile << EOF
${DIST} ${DELTA} ${NPT} 0.0 0.0
EOF
sprep96  -M ${VMOD} -HS 0 -HR 0 -L -R -NMOD ${NMODE} -d ddfile
sdisp96
sregn96 -NOQ
slegn96 -NOQ
FNYQ=`echo $DELTA | awk '{print 0.5/$1}' `

#####
#    make plot of the dispersion of the form
#    LC RC
#    LU RU
#####
rm -fr S?EGN?.PLT
rm -f DISP.PLT
sdpegn96 -L -C  -XLIN -YLIN -X0 2.0 -Y0 8 -XLEN 6 -YLEN 6 -YMIN 0 -YMAX ${VMAX} -XMIN 0.0 -XMAX ${FNYQ}
sdpegn96 -L -U  -XLIN -YLIN -X0 2.0 -Y0 1 -XLEN 6 -YLEN 6 -YMIN 0 -YMAX ${VMAX} -XMIN 0.0 -XMAX ${FNYQ}
sdpegn96 -R -C  -XLIN -YLIN -X0 9.5 -Y0 8 -XLEN 6 -YLEN 6 -YMIN 0 -YMAX ${VMAX} -XMIN 0.0 -XMAX ${FNYQ}
sdpegn96 -R -U  -XLIN -YLIN -X0 9.5 -Y0 1 -XLEN 6 -YLEN 6 -YMIN 0 -YMAX ${VMAX} -XMIN 0.0 -XMAX ${FNYQ}
cat S?EGN?.PLT > DISP.PLT

#####
#    now make the synthetics
#    for each subsource
#        get source coordinates
#        get force orientation
#        make synthetic
#        use gsac to apply the force
#        open the synthetic using cut o 0 o TMAX
#        save
#    then stack the subsources
#####
count=1
while [ $count -lt ${NSRC} ]
do
        SRC=`echo $count | awk '{printf "%6.6d", $1}' `
        getsrc
#       echo $EVLA $EVLO $XS $YS
        getforce
#       echo $FN $FE $FD
        #####
        #   Y = north
        #   X = east
        #####
        getdistaz $XS $YS $XR1 $YR1
        DIST1=$DIST
        AZ1=$AZ
        BAZ1=$BAZ
        getdistaz $XS $YS $XR2 $YR2
        DIST2=$DIST
        AZ2=$AZ
        BAZ2=$BAZ
        TSHIFT=`echo $SRC $NSRC $TMAX | awk '{WIN=$3/$2; print ($1 -1.) *  WIN}'`
        ######
        # check to see that DIST1 > DMIN and DIST2 > XMIN
        #####
        ANS=`echo $DIST1 $DIST2  $DMIN | awk '{ if ( $1 >= $3 && $2 >= $3 ) print "YES" ; else print "NO" }' `
        if [ $ANS = "YES" ]
        then
                makesyn $DIST1 $AZ1 $BAZ1 $XS $YS $XR1 $YR1
                mv T.Z 1.${SRC}.Z
                mv T.N 1.${SRC}.N
                mv T.E 1.${SRC}.E
                makesyn $DIST2 $AZ2 $BAZ2 $XS $YS $XR2 $YR2
                mv T.Z 2.${SRC}.Z
                mv T.N 2.${SRC}.N
                mv T.E 2.${SRC}.E
        count=`expr $count + 1 `
        fi
done


#####
#    make the final stack These are the 900 second long windows for noise
#####
gsac << EOF
cut o o ${TMAX}
r 1.??????.E
stack relative
ch kcmpnm E
w 1.E.stk

r 1.??????.N
stack relative
ch kcmpnm N
w 1.N.stk

r 1.??????.Z
stack relative
ch kcmpnm Z
w 1.Z.stk
q
EOF


gsac << EOF
cut o o ${TMAX}
r 2.??????.E
stack relative
ch kcmpnm E
w 2.E.stk

r 2.??????.N
stack relative
ch kcmpnm N
w 2.N.stk

r 2.??????.Z
stack relative
ch kcmpnm Z
w 2.Z.stk
q
EOF




    

DOSTKZNE

As set up, this script processes each of the 900 seconds of continuous 3-component noise at the two stations in 10 second segments. Each segment is whitened int he frequency domain and amplitude adjusted usg an AGC operator. Then the 10 second segments are cross-correlated. After all of the noise is processed, the cross-correlations and reversed cross-correlations are saved, and stacked to form the interstation Green's functions  Z.correv, N.correv and E.correv.   Because the stations were aligned in the EW direction, The Z.correv will have the Rayleigh wave, the E.correv will have the Rayleigh wave and the N.correv will have the  Love wave signal.

These were processed using the command

do_mft -G -IG ?.correv

The -G ensures that the selected dispersion files will be named Z.correv.dsp  for the ground velocities and Z.correc.phv for the phase velocities if the Z.correv is processed. The -IG flag
indicates that these are empirical Green's functions which permits the determination of phase velocities from the waveforms.

Recall that do_mft invokes the programs sacmft96 to process the waveform. Besides creating a data file of possible dispersion values and a figure, the scripts MFT96CMP and PHV96CMP are
created to permit the plot of theoretical dispersion on top of the output of sacmft96.

loveu.png
lovec.png
Fig. 2. Results of processing the file N.correv which will give the Love waves. The theoretical dispersion is plotted as red lines on top of the output of sacmft96.
The colors in the group velocity plot indicate the amplitude of the signal. The phase velocities are estimated from the largest amplitude in the ground velocity plot at a given period. The many phase velocity curves arise because of the ambiguity of multiples of 2 π radians in the phase.



raylu.png
raylc.png
Fig. 3. Results of processing the file Z.correv which will give the Rayleigh waves. The theoretical dispersion is plotted as red lines on top of the output of sacmft96. The colors in the group velocity plot indicate the amplitude of the signal. The phase velocities are estimated from the largest amplitude in the ground velocity plot at a given period. The many phase velocity curves arise because of the ambiguity of multiples of 2 π radians in the phase.




Last changed December 31, 2017