Difference between revisions of "Ambient vibration array"

From GeopsyWiki
Jump to navigation Jump to search
 
(77 intermediate revisions by 5 users not shown)
Line 3: Line 3:
 
== Common pre-requisites for array processing - getting ready ==
 
== Common pre-requisites for array processing - getting ready ==
  
=== Loading signals ===
+
For any array processing of ambient vibration data there are some
 +
common steps / pre-requisites regarding your signal recordings (selection
 +
of simultaneous recordings, coordinate insertion, signal grouping).
 +
Learn about all these options/preparatory steps in the following
 +
sub-chapters (linking to other parts of this wiki)
  
=== View time series ===
+
=== Loading and viewing signals (time series) ===
  
=== Group signals ===
+
Before processing your ambient vibration data, you need to [[Loading and viewing signals|load your signal files]]
 +
into geopsy and build a database. For checking the data quality of recordings and potential timing
 +
problems while recording you might want to use the [[Geopsy: Graphic|graphic viewer]] and the [[Geopsy: Chronogram|chronogram viewer]].
 +
It is often also a good idea to compute [[Spectral_amplitudes|Spectra]] and or [[H/V_spectral_ratio|H/V spectral ratios]] (in case of 3C data) from your
 +
recordings to check the spectral energy content, problematic installations and ''1-D-ness'' of your region.
  
=== Insert / edit station coordinates ===
+
links to [[Loading_and_viewing_signals|loading and viewing signals]]
  
== F-K Toolbox (conventional f-k) ==
+
=== Group signals ===
  
=== Opening f-k toolbox - overview ===
+
An array is defined by grouping a set of stations with common time base and recording simultaneously the properties
 +
of the wavefield at different, but closely co-located positions in space. After having built a database with many recordings, you may
 +
group individual arrays from the full set of signals. Both single component and three component arrays can be bundled
 +
together causing the array processing toolboxes to work on single or three component data, respectively.
  
[[Image:LepRing01OnlyVertical.png|thumb|right|400px|Signal viewer for one group of simultaneously recorded waveforms ready for array processing]]
+
For details on how to group signals for array processing take a look in [[Geopsy: Groups|groups]]
  
If you have followed all steps above you should have now a set of simultaneously recorded waveforms
+
=== Insert / edit station coordinates ===
building a group and containing coordinates for each waveform. You can view this group of waveforms in the
 
data viewer and should get a picture like the one shown to the right.
 
  
Now you can open the conventional frequency wavenumber toolbox by pushing the plugin icon [[Image:FKToolboxPluginIcon.png|32px]].
+
The final step in the preparation of your data set consists in attaching the stations coordinates to your
 +
recordings. Note that some [[Geopsy:_Supported_file_formats|waveform formats]] store the station coordinates and are eventually already included
 +
while importing the signal files. For most formats, however, this will be not the case. So, prepare you coordinate information
 +
from yor field book or your measurement devices (theodulite / (D)GPS device / etc.) and prepare some [[SciFigs: Multi-column Parser|multi-column ASCII file]]
 +
which can be read into [[geopsy]]. Details how to set coordinates can be found in section [[Geopsy: Set receivers|set receivers]]
  
The f-k toolbox should open and is now attached to the [[Geopsy: Signal Viewer|signal viewer]]. No other toolbox will now open with
 
the signals you have loaded to the signal viewer. In order to process the same data set by another
 
tool, you have to open them in a different signal viewer window. Note also, that once you close either
 
the f-k toolbox window or the signal viewer, there is no way to undo this action - you will have to start
 
again by opening a new signal viewer and attaching the f-k toolbox to it. Alternatively you can drag your
 
group directly into the f-k toolbox icon. Geopsy will automatically open a new exclusively attached signal viewer
 
for you. 
 
  
<br style="clear: both"/>
+
Since July 2018, high resolution and conventional FK toolboxes have been merged into a single FK toolbox. The text below should be reworked.
  
[[Image:FKToolboxTimeTab.png|thumb|left|300px|Conventional f-k toolbox]]
+
FK processing can be achieved inside Geopsy graphical interface or through a command line [[geopsy-fk]].
  
Let's now have a close look to the individual parameters to be set in the f-k toolbox. There are four
+
== F-K Toolbox (conventional f-k) ==
main groups of parameters to be considered.
 
  
* Pre-processing parameters for excluding data from processing.
+
The conventional frequency wavenumber technique as implemented in '''Geopsy''' is based
* Processing parameters for selecting time windows on the remaining data.
+
on the simple idea of [[Wikipedia:Sensor array#Delay-and-sum_beamforming|''delay and sum'' (or ''shift and sum'')]].
* Processing parameters for selecting narrow frequency bands for processing
+
This technique may be effectuated equivalently in time domain or frequency domain.
* Parameters related to the output of results from f-k processing.
+
In Geopsy we follow the frequency domain approach, as it is the most convenient and effective
 +
way to use this approach for determination of frequency dependent apparent velocity estimation
 +
(i.e dispersion curve estimation under the assumption of the wave field being composed of
 +
surface waves only).  
  
These parameters are laid out in different parts of the f-k toolbox widget.
+
The simultaneous waveform recordings of a group of spatially distributed stations are analyzed
Note, that the toolbox has three main tabs: ''Time'', ''Processing'' and ''Status''.  
+
in many narrow (mostly overlapping) frequency bands for individual analysis windows cut from the
 +
overall recordings. For each analysis window and frequency band, a grid search is performed
 +
in the wavenumber domain to effectively find the propagation properties of the most coherent
 +
and/or powerful plane wave arrival in the analysis window. Given the assumption of surface
 +
waves dominating the wave field, the apparent velocity equals the phase velocity of the
 +
surface wave at this particular frequency.
  
The pre-processing parameters for excluding data from processing follow the same argumentation as the
+
Details of the signal processing can be found in the [[Array signal processing|array signal processing]]  
data selection for other toolboxes like [[Geopsy: HV Toolbox|H/V processing]] and can be found in the
+
page of this wiki or following one of the links below in section [[#See also|see also]]
''Time tab''.
 
  
E.g. one can apply the so-called [[Geopsy: Anti Trigger|anti-trigger]]
+
For learning about the detailed use of the f-k toolbox (geopsy plugin) in a tutorial like fashion, please [[FK|follow this link to FK]].
approach for removing transient signals from the continuous recordings (i.e. excluding those time chunks from
 
processing). The anti-trigger approach can be applied to raw or filtered data. For the detailed
 
settings please look into the description of [[Geopsy: Anti Trigger|anti-trigger]].
 
Another, second, approach for excluding data from processing is provided via the [[Geopsy: Bad sample|bad sample]] approach.
 
For details in setting parameters with the bad sample approach check [[Geopsy: Bad sample|this page]].
 
  
<br style="clear: both"/>
+
== High resolution frequency wavenumber Toolbox (Capon's method) ==
 
 
=== Parameter settings ===
 
 
 
==== Parameters for time window settings ====
 
 
 
[[Image:FKTimeLimits.png|thumb|left|250px|Time limits layout group]]
 
There are two distinct inputs in the f-k toolbox for selecting time windows for processing.
 
First, it is necessary to specify start and end times of the overall data window for processing
 
by selecting entries ''From'' and ''To'' in the ''Time limits'' part of the ''Time tab'' of the fk-toolbox.
 
 
 
<br style="clear: both"/> 
 
 
 
There are three options each for the
 
drop boxes ''From'' and ''To''. For the ''From'' entry you can select among [[Image:FKFromDropList.png|100px]].
 
 
 
* this time (entries in form xx '''h''' xx '''m''' xx.xxx '''s''')
 
* T0 (selects earliest common start time of all traces in group)
 
* during (selects time duration before end of the signals - either specified with ''End'' or ''this time'' as explained below)
 
 
 
For the ''To'' entry there are the options [[Image:FKToDropList.png|100px]]:
 
 
 
* this time (entries in form xx '''h''' xx '''m''' xx.xxx '''s''')
 
* End (selects latest common end time of all traces in group)
 
* during (selects time duration from start time as given above - either specified by ''this time'' or ''T0'')
 
 
 
Note: a selection of time limits as ''From: during'' && ''To: during'' is NOT valid.
 
 
 
[[Image:FKGeneralTab.png|thumb|left|300px|General tab]]
 
After selecting the overall data window for processing you finally have to specify the window length for
 
each single analysis window to be processed. You can find this ''Length'' parameter in the ''General'' tab inside the
 
''Time tab''. There are again three options provided in a drop box [[Image:FKWindowLength.png|60px]].
 
 
 
* Exactly (takes exactly the time window length as provided in the corresponding editable box - unit is in seconds).
 
* At least (takes at least the time window length as given in the editable box, if possible it takes more - unit is in seconds).
 
* Freq. dep. (read: frequency dependent. Scales window length to the center frequency of each frequency band to be processed. The number given in the editable box corresponds to a window length in number of cycles for each frequency band. ''T'' is the center period of the band).
 
 
 
Note: the most common setting for the ''Length'' parameter is ''Freq. dep.'' with values between 20.00 to 50.00 meaning that 20 to 50 cycles
 
of the central siganl period per frequency bands are used as window length. Example: for a frequency band with center frequency 2.5 Hz, the central period is 1/2.5 = 0.4 s. Therefore, with a frequency dependent ''Length'' setting of ''30.00 T'', the window length will be chosen as 12.00 s.
 
  
<br style="clear: both"/>
+
The high resolution frequency wavenumber (hrfk) algorithm can be viewed as a '''''generalized''''' beamforming algorithm. It is implemented in '''Geopsy''' following the ideas of Capon (1969)
 +
<ref name="Capon (1969)"> Capon, J., High-Resolution Frequency-Wavenumber Spectrum Analysis, Proceedings of the IEEE, 57, No. 8, 1408-1419, 1969.</ref>,
 +
using an auto-adaptive (optimal) complex spatial weighting scheme for analysis of narrowband stationary signals. It is one of the most common and preferred frequency wavenumber techniques applied to
 +
ambient vibration analysis. 
  
==== Parameters for frequency band selection ====
+
The Capon beamformer relies on a very simple formula:
  
[[Image:FKToolboxProcessingTab.png|thumb|left|300px|Processing tab of f-k toolbox]]
 
The selection of the frequency bands to be processed is driven by parameter inputs inside the ''Processing'' tab of the f-k toolbox.
 
In the ''Frequency sampling'' layout group you can specify the minimum and the maximum (central) frequencies for processing.
 
Between these two limits, the frequency axis will then be sampled according to the settings in the next line [[Image:FKFreqSampling.png|250px]].
 
 
There are two options to be specified, namely ''Step'' and ''Number of Samples''. For the ''Step'' parameter, there are two options to be selected
 
from a drop box: ''Log'' or ''Linear''. You might have guessed the meaning right away. When choosing ''Linear'' then the frequency axis is sampled
 
linearly between the given limits and exactly ''Number of Samples'' will be distributed along the axis. Choosing ''Log'' accordingly samples the frequency axis logarithmically, i.e. frequency bands are ''close'' spaced at lower frequencies.
 
 
Note: the sampled frequencies are taken as '''center frequencies''' in the processing scheme. Therefore, it is necessary to specify additionally a bandwidth
 
so that the processing can be performed in a finite but narrow frequency band. The bandwidth parameter is chosen as relative half-bandwidth in the ''FK-gridding'' layout group in the ''processing tab'' [[Image:FKGridding.png|250px]].
 
 
The frequency limits for each frequency band are computed as
 
 
<math>
 
<math>
[(1-bw)\cdot f_c,(1+bw)\cdot f_c]
+
BP_{Capon}(\omega,\vec{k}) = \frac{1}{\vec{e}^\dagger(\omega,\vec{k})\underline{R}^{-1}(\omega)\vec{e}(\omega,\vec{k})}
 
</math>
 
</math>
where <math>f_c</math> is the central frequency and <math>bw</math> is the bandwidth parameter.
 
  
Example: a setting of 0.10 for the bandwidth, a ''Linear'' spacing of 10 samples along the frequency axis from 1 Hz to 10 Hz will create the following frequency bands:
+
where <math>\underline{R}^{-1}(\omega)</math> is simply the inverse of the cross spectral matrix estimate and <math>\vec{e}</math>
{| border="1" cellpadding="2" cellspacing="0"
+
is the so-called steering vector summarizing the shift times for the harmonic plane wave with wavenumber vector <math>\vec{k}</math>
|+ align="bottom" style="color:#e76700;" |Frequency band creation using ''Step'': Linear, ''Bandwidth'': 0.10, ''From'': 1 Hz, ''to'': 10 Hz
+
for each station within the array.The <math>\dagger</math> symbol presents the conjugate transpose operation. In the same notation,
! Freq. band index
+
the conventional f-k algorithm would be formulated as:
! Lower frequency limit
 
! Center frequency
 
! Upper frequency limit
 
|-
 
| 1
 
| 0.9 Hz
 
| 1.0 Hz
 
| 1.1 Hz
 
|-
 
| 2
 
| 1.8 Hz
 
| 2.0 Hz
 
| 2.2 Hz
 
|-
 
| 3
 
| 2.7 Hz
 
| 3.0 Hz
 
| 3.3 Hz
 
|-
 
| 4
 
| 3.6 Hz
 
| 4.0 Hz
 
| 4.4 Hz
 
|-
 
| 5
 
| 4.5 Hz
 
| 5.0 Hz
 
| 5.5 Hz
 
|-
 
| 6
 
| 5.4 Hz
 
| 6.0 Hz
 
| 6.6 Hz
 
|-
 
| 7
 
| 6.3 Hz
 
| 7.0 Hz
 
| 7.7 Hz
 
|-
 
| 8
 
| 7.2 Hz
 
| 8.0 Hz
 
| 8.8 Hz
 
|-
 
| 9
 
| 8.1 Hz
 
| 9.0 Hz
 
| 9.9 Hz
 
|-
 
| 10
 
| 9.0 Hz
 
| 10.0 Hz
 
| 11.0 Hz
 
|}
 
  
==== Parameters influencing the sampling of the wavenumber space ====
 
 
The basic idea of f-k processing consists of delaying the observed recordings at different stations according
 
to a particular horizontal wavenumber vector <math>\vec{k}=(k_x,k_y)^T</math> and computing the coherence and/or beam power of the
 
shifted stacked output of all array stations.
 
 
* Beampower:
 
<math>
 
BP(\omega,\vec{k}) = \sum_{\omega=\omega_l}^{\omega=\omega_h} \left[\sum_{i=1}^{i=N} X_i(\omega)\exp(j\vec{k}\vec{r}_i)\right]^2
 
</math>
 
* Semblance:
 
 
<math>
 
<math>
Semb(\omega,\vec{k}) = \frac{\sum_{\omega=\omega_l}^{\omega=\omega_h} \left[\sum_{i=1}^{i=N} X_i(\omega)\exp(j\vec{k}\vec{r}_i)\right]^2}
+
BP_{conv.}(\omega,\vec{k}) = \vec{e}^\dagger(\omega,\vec{k})\underline{R}(\omega)\vec{e}(\omega,\vec{k})
{N \sum_{\omega=\omega_l}^{\omega=\omega_h}\sum_{i=1}^{i=N} \left[X_i(\omega)\exp(j\vec{k}\vec{r}_i)\right]^2}
 
 
</math>
 
</math>
 
<math>X_i(\omega)</math> is the observed record at station <math>i</math> at frequency <math>\omega</math>. The station coordinates
 
in the array plane (usually earth's surface) are given by <math>\vec{r}_i</math>.
 
By testing many wavenumber vectors in the wavenumber plane one tries to find those wavenumbers which maximize the
 
array output. The wavenumber found corresponds to a plane wave along the surface crossing the seismometer array.
 
So, the problem at hand is a general optimization problem trying to find the maximum of a function of two parameters
 
(here <math>(k_x,k_y)^T</math>). In order to optimize the grid search approach, the f-k toolbox relies on the following sampling strategy.
 
 
 
<br style="clear: both"/>
 
 
=== Saving results / contents of output file ===
 
 
[[Image:FKOutput.png|thumb|right|250px|Power spectrum maxima group]]
 
 
There will be two output files created. The name of the output file and the location (path) for storing the output files is specified in the
 
lowest layout group of the ''Processing'' tab titled ''Power spectrum maxima'' with the option ''Output file''. Pushing the button to the right of the editable textbox entry will open a file browser and allows you to navigate to some location in your directory tree and specify an output file name.
 
 
The extension of the main output file is '''.max'''. Additionally an other output file is created containing the processing parameters and information created during
 
the processing run. The extension of this file is set to '''.log'''. Parameters can be reloaded from an existing ''.log'' file in the ''Load parameters'' option form the ''Time tab'' of the f-k toolbox [[Image:FKLoadParams.png|250px]].
 
 
<br style="clear: both"/>
 
 
# File generated by Geopsy, FK processing
 
# The process log is saved in file lep_ring01.log
 
# Header
 
# GRID_LAYOUT 0
 
# GRID_TYPE 0
 
# NPHI 72
 
# GRID_RESOL 250
 
# GRID_MAX 10
 
# NUM_BANDS 100
 
# Number of freq bands: 100
 
# Band 0 lower 0.45 center 0.5 upper 0.55
 
# Band 1 lower 0.468138 center 0.520153 upper 0.572169
 
# Band 2 lower 0.487007 center 0.541119 upper 0.595231
 
# Band 3 lower 0.506636 center 0.562929 upper 0.619222
 
# Band 4 lower 0.527057 center 0.585619 upper 0.644181
 
# Band 5 lower 0.548301 center 0.609223 upper 0.670146
 
# Band 6 lower 0.570401 center 0.633779 upper 0.697157
 
...............
 
...............
 
# Band 95 lower 19.2104 center 21.3449 upper 23.4794
 
# Band 96 lower 19.9847 center 22.2053 upper 24.4258
 
# Band 97 lower 20.7903 center 23.1003 upper 25.4103
 
# Band 98 lower 21.6282 center 24.0314 upper 26.4345
 
# Band 99 lower 22.5 center 25 upper 27.5
 
# seconds from start | cfreq | slow | az | math-phi | semblance | beampow
 
30378.9 0.940916 5.15376 21.9177 68.0823 0.491473 63.9522
 
30442.7 0.940916 3.18228 38.3304 51.6696 0.210185 65.6729
 
30474.6 0.940916 3.85247 121.714 328.286 0.298423 69.8839
 
30121.3 0.978841 4.66624 276.785 173.215 0.22663 70.1559
 
30366.6 0.978841 6.27949 18.5945 71.4055 0.188424 71.1481
 
30397.3 0.978841 6.79274 263.773 186.227 0.177627 63.9104
 
30354.7 1.0183 5.79576 65.0153 24.9847 0.192576 71.5275
 
30384.2 1.0183 6.28673 265.305 184.695 0.604842 71.942
 
30458.6 0.978841 4.57822 4.44118 85.5588 0.23338 69.5334
 
...........
 
...........
 
 
The contents of the output file are best explained by having a look into the plain ascii file. Header lines are marked with a ''#'' as the first character of the line. The header entries of the max files provide the list of frequency bands that were processed. Other entaa,ries are obsolete but are kept for compatibility with former file formats used in other codes during the [http://sesame-fp5.obs.ujf-grenoble.fr/index.htm SESAME project].
 
 
The more interesting lines of the the '''.max''' output file start right below the header lines. The final header line describes the entries in the different columns:
 
 
* Time from start of processed data window in seconds
 
* Center frequency of the processed data window
 
* Estimated Slowness in s/km (this unit is also kept for compatibility reasons)
 
* Azimuth of propagation direction (direction of wavenumber/slowness vector) as measured from North via East
 
* Same angle as before, but now given in mathematical orientation (measured from ''East'' via ''North'')
 
* Semblance coefficient
 
* Beampower value given in dB ( <math>10\cdot \log(BP(\omega,\vec{k}))</math>
 
 
Log file - different sections
 
 
section init
 
 
### Init Log ###
 
*********** vertical ***********
 
Add signal id 3 to component Vertical of station WA_WAU01 at 0.940538 -0.0688223 43.6439
 
Add signal id 6 to component Vertical of station WA_WAU02 at 5.32643 4.2312 43.7266
 
Add signal id 9 to component Vertical of station WA_WAU03 at 5.72042 -4.39827 42.5163
 
Add signal id 12 to component Vertical of station WA_WAU04 at 0.230693 -8.01841 42.9798
 
Add signal id 15 to component Vertical of station WA_WAU05 at -2.61856 -4.15918 43.5125
 
Add signal id 18 to component Vertical of station WA_WAU06 at -6.41618 0.965644 43.5362
 
Add signal id 21 to component Vertical of station WA_WAU07 at -3.53215 4.06034 43.6112
 
Add signal id 24 to component Vertical of station WA_WAU08 at 0.34882 7.3875 43.6973
 
Found 8 different stations
 
### End Init Log ###
 
 
Parameter section - note this can be loaded via the ''Load Params'' button.
 
 
### Parameters ###
 
WINDOW MIN LENGTH (s) = 50
 
WINDOW MAX LENGTH (s) = 30
 
WINDOW LENGTH TYPE (at least/exactly/freq. dep.) = freq. dep.
 
DO BAD SAMPLE TOLERANCE (y/n) = n
 
BAD SAMPLE TOLERANCE (s) = 0
 
DO WINDOW OVERLAP (y/n) = n
 
WINDOW OVERLAP (%) = 5
 
DO BAD SAMPLE THRESHOLD (y/n) = n
 
BAD SAMPLE THRESHOLD (%) = 99
 
ANTI-TRIGGERING ON RAW SIGNAL (y/n) = n
 
USED RAW COMPONENTS = y, n, y, y, y, y, y, y, y, y
 
RAW STA (s) = 1
 
RAW LTA (s) = 30
 
RAW MIN SLTA = 0.2
 
RAW MAX SLTA = 2.5
 
MINIMUM FREQUENCY = 0.5
 
MAXIMUM FREQUENCY = 25
 
INVERSED FREQUENCY = n
 
SAMPLES NUMBER FREQUENCY = 100
 
SAMPLING TYPE FREQUENCY (0=log, 1=linear)= 0
 
FROM TIME TYPE = 0
 
FROM TIME TEXT = 8h21m0s
 
TO TIME TYPE = 0
 
TO TIME TEXT = 9h29m0s
 
MIN K = 0.015
 
MAX K = 0.5
 
MIN V = 100
 
FREQ BAND WIDTH = 0.1
 
N MAXIMA = 3
 
OUTPUT FILE = /home/mao/GEOPSY_DOC_WORKSHOP/lep_ring01.max
 
### End Parameters ###
 
 
Process Log section
 
 
### Process Log ###
 
Process started at 2010-03-08 10:54:47
 
Adding window from 30060 to 30120 s.
 
Adding window from 30120 to 30180 s.
 
Adding window from 30180 to 30240 s.
 
Adding window from 30240 to 30300 s.
 
Adding window from 30300 to 30360 s.
 
Adding window from 30360 to 30420 s.
 
Adding window from 30420.1 to 30480.1 s.
 
Adding window from 30480.1 to 30540.1 s.
 
Adding window from 30540.1 to 30600.1 s.
 
Adding window from 30600.1 to 30660.1 s.
 
Adding window from 30660.1 to 30720.1 s.
 
Adding window from 30720.1 to 30780.1 s.
 
Adding window from 30780.1 to 30840.1 s.
 
Adding window from 30840.1 to 30900.1 s.
 
Frequency 1/100 0.5
 
Min Window length 60 seconds
 
Max Window length 60 seconds
 
  14 Time windows
 
.......
 
Adding window from 30919.1 to 30920.3 s.
 
Adding window from 30920.3 to 30921.5 s.
 
Adding window from 30921.5 to 30922.7 s.
 
Frequency 100/100 25
 
Min Window length 1.2 seconds
 
Max Window length 1.2 seconds
 
  713 Time windows
 
Process run in 00:01:20
 
Process ended at 2010-03-08 10:56:08
 
### End Process Log ###
 
 
=== Graphical display of f-k results using max2curve ===
 
 
[[max2curve]]
 
 
== High resolution frequency wavenumber Toolbox (Capon's method) ==
 
 
The high resolution frequency wavenumber (hrfk) algorithm follows the ideas of Capon (1969)
 
<ref name="Capon (1969)"> Capon, 1969. xxx ,IEEE Proceedings, , , .</ref>
 
 
. The corresponding toolbox that can be called using the following
 
plugin icon [[Image:HRFKPluginIcon.png|32px]] can be used straightforward in nearly the same way
 
as the conventional f-k toolbox.
 
 
There is just one more (optional!) parameter that can be set. The calcuation for the hrfk algorithm involves the inversion of
 
the cross spectral matrix estimate from the data. In order to guarantee a numerically stable solution for the Matrix inversion,
 
a regularization parameter (damping constant) can be provided to load the diagonal of the cross spectral matrix before inversion.
 
 
Note: by doing so (using a damping constant), you will smooth out your wavenumber power spectra and eventually not obtain
 
much better results compared to the conventional beamforming approach.
 
 
 
   
 
   
 +
The required computation of the inverse of the cross-spectral matrix results in some potential instabilities for the HRFK estimates.
  
=== Saving results / contents of output file ===
+
The corresponding toolbox can be opened using the icon [[Image:HRFKPluginIcon.png|32px]]. It is nearly equivalent to the [[FK|conventional frequency wavenumber]]
 +
toolbox and the processing flow is exactly equivalent to the conventional technique.
  
=== Graphical display of f-k results using max2curve ===
+
For learning about the detailed use of the HRFK toolbox (geopsy plugin) in a tutorial like fashion, please [[HRFK| follow this link to HRFK]]
  
 
== Modified Spatial Autocorrelation (MSPAC) Toolbox ==
 
== Modified Spatial Autocorrelation (MSPAC) Toolbox ==
  
[[Image:LepRing01OnlyVertical_SPAC_VIEWALL_NEW.png|thumb|right|600px]]
+
The Modified Spatial Autocorrelation (MSPAC) was introduced by Bettig et al. (2001) <ref name="Bettig et al. (2001)"> Bettig B., P.-Y. Bard, F. Scherbaum, J. Riepl, F. Cotton, C. Cornou, D. Hatzfeld, 2001. Analysis of dense array measurements using the modified spatial auto-correlation method (SPAC). Application to Grenoble area., Boletin de Geofisica Teoria e Applicata, 42, 3-4, 281-304.</ref> after pioneer paper of Aki (1957)<ref name="Aki (1957)."> Aki, K., 1957. Space and Time Spectra of Stationary Stochastic Waves, with Special Reference to Microtremors, Bull. Earthq. Res. Inst. Tokyo, 35, 415-457. </ref> and allows to compute average spatial autocorrelation coefficients for any arbitary array configurations. As SPAC technique, MSPAC relies on a stochastic ambient noise wavefield stationary in both time and space. Aki (1957) showed, that, given this assumption, the existing relation between the spectrum densities in space and time can be used to derive the following expression for a plane wave narrow-band filtered around <math>\omega_0</math>:
The Modified Spatial Autocorrelation (MSPAC) was introduced by Bettig et al. (2001) <ref name="Bettig et al. (2001)"> Bettig B., P.-Y. Bard, F. Scherbaum, J. Riepl, F. Cotton, C. Cornou, D. Hatzfeld, 2001. Analysis of dense array measurements using the modified spatial auto-correlation method (SPAC). Application to Grenoble area., Boletin de Geofisica Teoria e Applicata, 42, 3-4, 281-304.</ref> after pioneer paper of Aki (1957)<ref name="Aki (1957)."> Aki, K., 1957. Space and Time Spectra of Stationary Stochastic Waves, with Special Reference to Microtremors, Bull. Earthq. Res. Inst. Tokyo, 35, 415-457. </ref>.  This method allows computing spatial autocorrelation coefficients for any arbitary array layouts.
 
  
After viewing the group of waveforms in the data viewer, open the MSPAC toolbox by pushing the following plugin icon  [[Image:SPACToolboxPluginIcon.png|18px]]. Alternatively you can drag your group directly into the SPAC toolbox icon. The SPAC toolbox should now be attached to the signal viewer.
+
<math>
Four windows should then opened as displayed in the figure on the right: (1) a signal time viewer, (2) a map displaying array station locations, (3) a map displaying co-array stations location and (4) the SPAC processing toolbox. This processing toolbox is composed of four sub-menus: "Rings" sub-menu is used to design rings from the co-array; "Time" is used to define begin and end time to be processed, time window length, anti-triggering options as described in section '''XXX'''; "Output" menu is used to set up the output frequency band together with the sampling scale type (linear, log) and number of frequencies, the output filename; "Status" menu provides information about computing status. All these sub-menus are described hereafter.
+
\overline{\rho(\omega_0 , r)} = \frac{1}{\pi}\int^{0}_{\pi}\rho(\omega_0,r,\varphi)d(\varphi)= J_0(\frac{\omega r}{c(\omega_0)})
 +
</math>
  
<br style="clear: both"/>
+
<math>\overline{\rho(\omega_0 , r)}</math> represents the averaging over azimuth of spatial autocorrelations  <math> {\rho(\omega_0 , r, \varphi)} = cos(\frac{\omega_0 r}{c(\omega_0)}cos(\theta-\varphi)) </math> where <math>{\theta}</math> is the wave azimuth and <math>{\varphi}</math> the direction azimuth between stations pairs. 
  
[[Image:LepRing01OnlyVertical_SPAC_SELECT_RINGS_FIRST_STEP.png|thumb|left|500px|Steps involved in the design of the rings.]]
 
  
=== Defining rings ===
 
  
In "rings" sub-menu:
+
Application of the SPAC technique (as well as further derived techniques like ESAC) requires perfect shaped arrays (circular, semi-circular, nested triangles) which may be difficult to achieve in urban environment. To overcome these difficulties, Bettig et al. (2001) <ref name="Bettig et al. (2001)"> Bettig B., P.-Y. Bard, F. Scherbaum, J. Riepl, F. Cotton, C. Cornou, D. Hatzfeld, 2001. Analysis of dense array measurements using the modified spatial auto-correlation method (SPAC). Application to Grenoble area., Boletin de Geofisica Teoria e Applicata, 42, 3-4, 281-304.</ref> suggested to use the co-array - which is defined as the set of all possible combinations of two array sensors (Haubrich, 1968 <ref name="Haubrich (1968)"> Haubrich, R.A., 1968. Array Design, Bull. seism. Soc. Am., 58(3), 977–991. </ref>) - to divide the array into several semicircular sub-arrays. Each sub-array (hereafter called ring) is thus composed of several sensors pairs. To account for ring thickness, an azimuthal and radial integration is then needed to compute averaged spatial autocorrelation values <ref name="Bettig et al. (2001)"> Bettig B., P.-Y. Bard, F. Scherbaum, J. Riepl, F. Cotton, C. Cornou, D. Hatzfeld, 2001. Analysis of dense array measurements using the modified spatial auto-correlation method (SPAC). Application to Grenoble area., Boletin de Geofisica Teoria e Applicata, 42, 3-4, 281-304.</ref>:
* Push on [[Image:SPACToolbox_AddIcon.png|50px]] button  to add a ring
 
* Specify the inner and outer circles radii of the ring (see step (1) in the figure on the left)
 
* Press on [[Image:SPACToolbox_OptimizeIcon.png|50px]] buttom to adjust the inner and outer circles radii in order to fit at best sensors location (see step (2) in the figure on the left). This ring is composed of 10 pairs of sensors as specified in the "Pairs" column.
 
* Define the next rings (see step (3) in the figure on the left). Note that you can associate a specific color to each ring.
 
* Note that you add or remove specific ring by pushing on the [[Image:SPACToolbox_AddIcon.png|50px]] or [[Image:SPACToolbox_RemoveIcon.png|50px]] buttons
 
* Once done, click on [[Image:SPACToolbox_SaveIcon.png|50px]] to save your rings design.  
 
[[Image:LepRing01OnlyVertical_SPAC_SAVE_RINGS.png|thumb|center|300px]]
 
* The rings design can next be loaded by pushing on the [[Image:SPACToolbox_LoadIcon.png|50px]] button
 
* the rings design are saved in ascii format with the following format: inner circle radius (in meters), outer circle radius (in meters), RGB color code associated to the ring
 
  
# MinRadius    MaxRadius  Red  Green  Blue
+
<math>
  4.23 ...
+
\overline{\rho_{r_1,r_2}(\omega)} = \frac{2}{r_2^2-r_1^2} \int^{r_1}_{r_2}r.J_0(\frac{\omega r}{c(\omega)})dr = \frac{2}{r_2^2-r_1^2}\frac{c(\omega)}{\omega}[r.J_1(\frac{\omega r}{c(\omega)})]^{r_1}_{r_2}
 
+
</math>
[[Image:LepRing01OnlyVertical_SPAC_EDIT_RINGS.png|thumb|center|400px|Rings definition format]]
 
 
 
 
 
<br style="clear: both"/>
 
 
 
=== Defining time limits ===
 
 
 
In "Time" sub-menu:
 
[[Image:SPACToolbox_TimeLimits.png|thumb|right|300px]]
 
* Select the time limits to be processed. You can either use the begin '''To''' and end time '''End''' of the time series, either specify the begin and end time ''' This Time ''' or specify the time duration '''Duration'''
 
 
 
<br style="clear: both"/>
 
 
 
[[Image:SPACToolbox_WindowLength.png|thumb|right|300px]]
 
* Specify the window length that will be used for processing. You can either choose a window length that is inversely proportionnal to the frequency of interest '''Freq. Dep.''', a fixed time window length '''Exactly''' or a variable length limited by a minimum and a maximumm '''At Least'''.
 
* By pushing on [[Image:SPACToolbox_Parameters.png|80px]] , you can also load the .log file containing all the input parameters saved from a previous computation (see '''Editing spatial autocorrelation curves''' for complementary information).
 
 
 
<br style="clear: both"/>
 
 
 
=== Defining the outputs and starting the computation ===
 
 
 
In "Output" sub-menu:
 
[[Image:SPACToolbox_FrequencyLimits.png|thumb|left|300px]]
 
* Choose the frequency band used for the computation
 
* Specify the number of frequencies and the sampling scale (constant on a linear or logarithmic scale).
 
<br style="clear: both"/>
 
[[Image:SPACToolbox_OutputName.png|thumb|left|300px]]
 
* Click on '''Dinver target''' to save your results in a .target file. Set the output filename.
 
* Then, click on [[Image:SPACToolbox_StartIcon.png|50px]] to start the computation. Status of the computation can be checked in the "Status" sub-menu as displayed below.
 
[[Image:SPACToolbox_ComputationStatus.png|thumb|center|200px]]
 
 
 
<br style="clear: both"/>
 
 
 
=== Graphical display of spatial autocorrelation curves / contents of output files ===
 
 
 
[[Image:SPACToolbox_LogFile.png|thumb|right|400px|Contents of the .log file after MSPAC computation.]]
 
The spatial autocorrelation coefficients and the parameters used for the computation are saved in a .target binary file and a .log ascii file, respectively. The .log file is saved in the same folder as the .target file. The figure on the right displays the content of the .log file.  The .target file is a binary file.
 
  
Graphical display of spatial autororrelation curves can be done by importing the .target file in the "target" menu of "Dinver" (CROSS-REFERENCE) or by using an utility called '''spac2disp''' that can be launched either from the windows Geopsy menu either from a command line
+
where <math> r_1 </math> and <math> r_2 </math> are the inner and outer radius of the ring, respectively.
  
[[Image:SPACToolbox_LaunchSpac2disp.png|thumb|center|400px]]
+
This method allows computing spatial autocorrelation coefficients for any arbitary array configurations. The design of rings results from a compromise between a number of sensors pair per ring as large as possible and a ring thickness as small as possible.
  
  
<br style="clear: both"/>
+
For learning the usage of the MSPAC Toolbox in detail in a tutorial like fashion, please [[MSPAC| follow this link to MSPAC]].
  
 +
== See also ==
  
 +
* [[Wikipedia:Array_signal_processing|Wikipedia page on array signal processing]]
 +
* [http://www.geopsy.org geopsy.org]
 +
* etc.
  
=== References ===
+
== References ==
 
<references/>
 
<references/>

Latest revision as of 07:37, 10 July 2018


Common pre-requisites for array processing - getting ready

For any array processing of ambient vibration data there are some common steps / pre-requisites regarding your signal recordings (selection of simultaneous recordings, coordinate insertion, signal grouping). Learn about all these options/preparatory steps in the following sub-chapters (linking to other parts of this wiki)

Loading and viewing signals (time series)

Before processing your ambient vibration data, you need to load your signal files into geopsy and build a database. For checking the data quality of recordings and potential timing problems while recording you might want to use the graphic viewer and the chronogram viewer. It is often also a good idea to compute Spectra and or H/V spectral ratios (in case of 3C data) from your recordings to check the spectral energy content, problematic installations and 1-D-ness of your region.

links to loading and viewing signals

Group signals

An array is defined by grouping a set of stations with common time base and recording simultaneously the properties of the wavefield at different, but closely co-located positions in space. After having built a database with many recordings, you may group individual arrays from the full set of signals. Both single component and three component arrays can be bundled together causing the array processing toolboxes to work on single or three component data, respectively.

For details on how to group signals for array processing take a look in groups

Insert / edit station coordinates

The final step in the preparation of your data set consists in attaching the stations coordinates to your recordings. Note that some waveform formats store the station coordinates and are eventually already included while importing the signal files. For most formats, however, this will be not the case. So, prepare you coordinate information from yor field book or your measurement devices (theodulite / (D)GPS device / etc.) and prepare some multi-column ASCII file which can be read into geopsy. Details how to set coordinates can be found in section set receivers


Since July 2018, high resolution and conventional FK toolboxes have been merged into a single FK toolbox. The text below should be reworked.

FK processing can be achieved inside Geopsy graphical interface or through a command line geopsy-fk.

F-K Toolbox (conventional f-k)

The conventional frequency wavenumber technique as implemented in Geopsy is based on the simple idea of delay and sum (or shift and sum). This technique may be effectuated equivalently in time domain or frequency domain. In Geopsy we follow the frequency domain approach, as it is the most convenient and effective way to use this approach for determination of frequency dependent apparent velocity estimation (i.e dispersion curve estimation under the assumption of the wave field being composed of surface waves only).

The simultaneous waveform recordings of a group of spatially distributed stations are analyzed in many narrow (mostly overlapping) frequency bands for individual analysis windows cut from the overall recordings. For each analysis window and frequency band, a grid search is performed in the wavenumber domain to effectively find the propagation properties of the most coherent and/or powerful plane wave arrival in the analysis window. Given the assumption of surface waves dominating the wave field, the apparent velocity equals the phase velocity of the surface wave at this particular frequency.

Details of the signal processing can be found in the array signal processing page of this wiki or following one of the links below in section see also

For learning about the detailed use of the f-k toolbox (geopsy plugin) in a tutorial like fashion, please follow this link to FK.

High resolution frequency wavenumber Toolbox (Capon's method)

The high resolution frequency wavenumber (hrfk) algorithm can be viewed as a generalized beamforming algorithm. It is implemented in Geopsy following the ideas of Capon (1969) [1], using an auto-adaptive (optimal) complex spatial weighting scheme for analysis of narrowband stationary signals. It is one of the most common and preferred frequency wavenumber techniques applied to ambient vibration analysis.

The Capon beamformer relies on a very simple formula:

where is simply the inverse of the cross spectral matrix estimate and is the so-called steering vector summarizing the shift times for the harmonic plane wave with wavenumber vector for each station within the array.The symbol presents the conjugate transpose operation. In the same notation, the conventional f-k algorithm would be formulated as:

The required computation of the inverse of the cross-spectral matrix results in some potential instabilities for the HRFK estimates.

The corresponding toolbox can be opened using the icon HRFKPluginIcon.png. It is nearly equivalent to the conventional frequency wavenumber toolbox and the processing flow is exactly equivalent to the conventional technique.

For learning about the detailed use of the HRFK toolbox (geopsy plugin) in a tutorial like fashion, please follow this link to HRFK

Modified Spatial Autocorrelation (MSPAC) Toolbox

The Modified Spatial Autocorrelation (MSPAC) was introduced by Bettig et al. (2001) [2] after pioneer paper of Aki (1957)[3] and allows to compute average spatial autocorrelation coefficients for any arbitary array configurations. As SPAC technique, MSPAC relies on a stochastic ambient noise wavefield stationary in both time and space. Aki (1957) showed, that, given this assumption, the existing relation between the spectrum densities in space and time can be used to derive the following expression for a plane wave narrow-band filtered around :

represents the averaging over azimuth of spatial autocorrelations where is the wave azimuth and the direction azimuth between stations pairs.


Application of the SPAC technique (as well as further derived techniques like ESAC) requires perfect shaped arrays (circular, semi-circular, nested triangles) which may be difficult to achieve in urban environment. To overcome these difficulties, Bettig et al. (2001) [2] suggested to use the co-array - which is defined as the set of all possible combinations of two array sensors (Haubrich, 1968 [4]) - to divide the array into several semicircular sub-arrays. Each sub-array (hereafter called ring) is thus composed of several sensors pairs. To account for ring thickness, an azimuthal and radial integration is then needed to compute averaged spatial autocorrelation values [2]:

where and are the inner and outer radius of the ring, respectively.

This method allows computing spatial autocorrelation coefficients for any arbitary array configurations. The design of rings results from a compromise between a number of sensors pair per ring as large as possible and a ring thickness as small as possible.


For learning the usage of the MSPAC Toolbox in detail in a tutorial like fashion, please follow this link to MSPAC.

See also

References

  1. Capon, J., High-Resolution Frequency-Wavenumber Spectrum Analysis, Proceedings of the IEEE, 57, No. 8, 1408-1419, 1969.
  2. 2.0 2.1 2.2 Bettig B., P.-Y. Bard, F. Scherbaum, J. Riepl, F. Cotton, C. Cornou, D. Hatzfeld, 2001. Analysis of dense array measurements using the modified spatial auto-correlation method (SPAC). Application to Grenoble area., Boletin de Geofisica Teoria e Applicata, 42, 3-4, 281-304.
  3. Aki, K., 1957. Space and Time Spectra of Stationary Stochastic Waves, with Special Reference to Microtremors, Bull. Earthq. Res. Inst. Tokyo, 35, 415-457.
  4. Haubrich, R.A., 1968. Array Design, Bull. seism. Soc. Am., 58(3), 977–991.