Difference between revisions of "Geopsy-fk"

From GeopsyWiki
Jump to navigation Jump to search
 
(119 intermediate revisions by the same user not shown)
Line 1: Line 1:
 
== Signal database preparation ==
 
== Signal database preparation ==
  
[[Image:LepRing01OnlyVertical.png|thumb|right|400px|Signal viewer for one group of simultaneously recorded waveforms ready for array processing]]
+
This tutorial is based on experimental data recorded in Mirandola (Italy) during [https://interpacific.geopsy.org InterPACIFIC] project. The details of the experimental setup are available in Garofalo et al. (2016) <ref name="Garofalo(2016)">Garofalo, F., Foti, S., Hollender, F., Bard, P., Cornou, C., Cox, B. R., Ohrnberger, M., Sicilia, D., Asten, M., Di Giulio, G., Forbriger, T., Guillier, B., Hayashi, K., Martin, A., Matsushima, S., Mercerat, D., Poggi, V., Yamanaka, H. (2016). InterPACIFIC project: comparison of invasive and non-invasive methods for seismic site characterization. Part I: intra-comparison of surface wave methods, ''Soil Dynamics and Earthquake Engineering'', '''82''', 222-240, [https://doi.org/10.1016/j.soildyn.2015.12.010 10.1016/j.soildyn.2015.12.010]</ref>.
  
This step is performed with Geopsy graphical interface. If you have followed all steps for the preparation of a [[Geopsy: Database|database]], [[Geopsy: Set receivers|loading coordinates]]  
+
Download signals [https://interpacific.geopsy.org/data/AVA_MIRANDOLA_SAC.zip AVA in SAC] for Mirandola from [https://interpacific.geopsy.org InterPACIFIC]. Decompress the ZIP file. FK processing can access data in two ways:
and [[Geopsy: Groups|grouping]], you should have now a set of simultaneously recorded waveforms associated in a group and containing coordinates for each waveform.
+
* through a [[#Database|database]]
 +
* [[#Direct file access|directly reading]] the signal files.
  
Note that from release 3.4.0, it is possible to load signal files directly from the command line without creating a database. The file headers must have the correct absolute time and the sensor coordinates.
+
=== Database ===
  
Alternatively, you can download [[Media:lep_ring01_3C.gpy|a ready-to-go database file]] and the
+
If you have time, you can build your own database.  Alternatively, you can skip this tedious step by downloading [[Media:Mirandola.gpy|Mirandola.gpy]]. A [[#Using an existing database|minor adjustment]] is still necessary.
corresponding waveform files at [[Media:LEP20100224_RING01.tgz‎|this location]]. Note that you have to re-locate your
 
waveform files as described in detail [[Geopsy: Open%2C Close and Save|in this page]]. Or without a graphical interface:
 
  
   $ geopsy-fk -db Lep_ring01_3C.gpy -fix-signal-paths
+
==== Creating a database ====
 +
 
 +
To build a database, start by importing the signals in ''geopsy'' graphical user interface and by creating a group for each array. There are two options to import signals located in various sub-directories:
 +
 
 +
  $ geopsy MIRANDOLA_SAC/*/*.sac
 +
 
 +
or by launching ''geopsy'' and importing signals inside the graphical interface with the menu '''File/Import signals/File pattern'''. Each sub-directory is a synchronous array. Create a [[Geopsy:_Groups|group]] for each one, eventually with option '''Default array groups'''. Set the coordinates of the sensors with the menu '''Edit/[[Geopsy: Set receivers|Set receivers]]'''. The coordinates are provided in the same archive as the signals (files ''.geom''). Finally save the database to a new .gpy file.
 +
 
 +
==== Using an existing database ====
 +
 
 +
The section applies only if you intend to use a database created on another machine. A database file contains absolute paths to the signal files. On your computer the original paths certainly do not exist. Paths can be fixed inside ''geopsy'' graphical interface but it is also possible to do it through the command line, helpful if you are running on a remote server.
 +
 
 +
   $ geopsy-fk -db Mirandola.gpy -fix-signal-paths
 
   ----WARNING--- Opening Database...----
 
   ----WARNING--- Opening Database...----
   File '/home/mao/GEOPSY_DOC_WORKSHOP/RING01_SHORT/WA.WAU01..HHE.D.2010.056.000' does not exist.
+
   File '/home/wathelem/Sites/Mirandola/MIRANDOLA_SAC/MIR_C_135_405/MIR_C_135_405_CN01_E.sac' does not exist.
 
   The path may have changed. Would you like to manually select its new location?
 
   The path may have changed. Would you like to manually select its new location?
 
     1. Yes <-- default
 
     1. Yes <-- default
Line 20: Line 31:
 
   ?
 
   ?
  
Hit enter or answer '1', 'y' or 'yes'.
+
Hit enter or answer '1', 'y', 'yes' or just press ENTER.
  
 
   Show this message again? [y]/n
 
   Show this message again? [y]/n
  
Hit enter or answer 'y' or 'yes'.
+
Hit enter or answer 'y', 'yes' or just press ENTER.
  
 
   ---- Opening Database... ----
 
   ---- Opening Database... ----
     Filter: Signal file (WA.WAU01..HHE.D.2010.056.000)
+
     Filter: Signal file (MIR_C_135_405_CN01_E.sac)
 
     Current directory: /tmp
 
     Current directory: /tmp
     Please select file '/home/mao/GEOPSY_DOC_WORKSHOP/RING01_SHORT/WA.WAU01..HHE.D.2010.056.000' in its new location.
+
     Please select file '/home/wathelem/Sites/Mirandola/MIRANDOLA_SAC/MIR_C_135_405/MIR_C_135_405_CN01_E.sac' in its new location.
 
   File to open:
 
   File to open:
  
 
Provide the absolute or the relative path to the requested file(s). For instance:
 
Provide the absolute or the relative path to the requested file(s). For instance:
  
   RING01_SHORT/WA.WAU01..HHE.D.2010.056.000
+
   MIRANDOLA_SAC/MIR_C_135_405/MIR_C_135_405_CN01_E.sac
  
Files having similar paths are automatically translated. At the end the database is saved with your local paths. You can check that all paths are correct by re-starting the same command: you should get no output message. You can also view the signals in a graphical user interface:
+
Files having similar paths are automatically translated. At the end the database is saved with your local paths. You can check that all paths are correct by re-starting the same command: you should get no output message. Under Windows, you must use '\' as separator when specifying the new path. It is probably more efficient to do it with the graphical interface under this platform:
  
   $ geopsy Lep_ring01_3C.gpy
+
   $ geopsy Mirandola.gpy
  
 +
The list of available groups (name, starting and ending time, duration):
  
This database contains two predefined groups: ''vertical'' and ''3C''. Please test if you can [[Loading and viewing signals|view]]
+
  $ geopsy-fk -db Mirandola.gpy -groups
the signals: use [[Geopsy: Signal drag%26drop|drag and drop]] functionality and check the coordinate settings in the
+
  geopsy-fk: C_5_15-3C: 20130829100530.000000 20130829110530.000000 (3600 sec.)
[[Geopsy: Table|table view]] or [[Geopsy: Map|map view]]. If you display the predefined group ''vertical''
+
  geopsy-fk: C_5_15-Z: 20130829100530.000000 20130829110530.000000 (3600 sec.)
in a [[Geopsy: Graphic|graphic viewer]], you should obtain a picture similar to the one displayed on the right.
+
  geopsy-fk: C_15_45-3C: 20130829113830.000000 20130829125330.960000 (4500.96 sec.)
 +
  geopsy-fk: C_15_45-Z: 20130829113830.000000 20130829125330.960000 (4500.96 sec.)
 +
  geopsy-fk: C_45_135-3C: 20130829134300.000000 20130829145600.960000 (4380.96 sec.)
 +
  geopsy-fk: C_45_135-Z: 20130829134300.000000 20130829145600.960000 (4380.96 sec.)
 +
  geopsy-fk: C_26_78-3C: 20130830063630.000000 20130830075631.130000 (4801.13 sec.)
 +
  geopsy-fk: C_26_78-Z: 20130830063630.000000 20130830075631.130000 (4801.13 sec.)
 +
  geopsy-fk: C_135_405-3C: 20130829155630.000000 20130829175500.000000 (7110 sec.)
 +
  geopsy-fk: C_135_405-Z: 20130829155630.000000 20130829175500.000000 (7110 sec.)
  
Another way to get the list of groups is:
+
=== Direct file access ===
  
   $ geopsy-fk -db Lep_ring01_3C.gpy -groups
+
This option works only if the coordinates are set in the signal file headers which is not the case for proposed example. The quotes are mandatory to avoid the expansion of file names by the shell.
   geopsy-fk: Default groups
+
 
   geopsy-fk: Default groups/All signals
+
   $ geopsy-fk "MIRANDOLA_SAC/MIR_C_135_405/*.sac"
   geopsy-fk: Default groups/Temporary signals
+
  [...]
   geopsy-fk: Default groups/All files
+
   geopsy-fk: Relative coordinates with reference at (0.00 0.00).
   geopsy-fk: Default groups/Temporary files
+
   geopsy-fk: Null coordinates
   geopsy-fk: Default groups/Permanent files
+
   geopsy-fk: Error initializing group 'All signals'
   geopsy-fk: Default groups/By names
+
 
   geopsy-fk: Default groups/By components
+
To check the coordinates:
   geopsy-fk: Default groups/By receivers
+
 
   geopsy-fk: Default groups/By sources
+
   $ geoosy-fk "MIRANDOLA_SAC/MIR_C_135_405/*.sac" -coord
  geopsy-fk: vertical
+
  0 0 0 CN01
  geopsy-fk: 3C
+
   0 0 0 CN02
 +
   0 0 0 CN03
 +
   [...]
 +
 
 +
The option ''-utm'' can be used to assign coordinates to the signals. The accepted format is relatively rigid: utm_zone x y station_name. The file provided by InterPACIFIC project were edited under Windows, that is, the lines are all ended by CR and LF. Under Linux, a single LF is expected. The following script removes CR and re-orders the columns to fit the expected format.
 +
 
 +
   cat MIRANDOLA_GEOM/MIR_C_135_405.geom | sed "s/\r$//" |  \
 +
   awk '{if($2!="" && NR>1) {print "---", $2, $3, $1}}' > MIRANDOLA_GEOM/MIR_C_135_405.coord
 +
 
 +
You can check that the coordinates are properly assigned to the signals.
 +
 
 +
   $ geopsy-fk "MIRANDOLA_SAC/MIR_C_135_405/*.sac" -utm MIRANDOLA_GEOM/MIR_C_135_405.coord -coord
  
 
== Single component high resolution FK ==
 
== Single component high resolution FK ==
 +
 +
[[Image:Mirandola-capon-first.png|thumb|right|500px|Dispersion curve obtained with the vertical component high resolution FK (default parameters).]]
  
 
The simplest way to invoke ''geopsy-fk'' is:
 
The simplest way to invoke ''geopsy-fk'' is:
  
   $ geopsy-fk -db ../Lep_ring01_3C.gpy -group vertical
+
   $ geopsy-fk -db Mirandola.gpy -group C_135_405-Z
 +
 
 +
or
 +
 
 +
  $ geopsy-fk "MIRANDOLA_SAC/MIR_C_135_405/*.sac" -utm MIRANDOLA_GEOM/MIR_C_135_405.coord -set PROCESS_TYPE=Capon
 +
 
 +
according to the method chosen in the previous section. In the second option, you must specify the process type to ''Capon'' (single component high resolution FK) unless you want a three component process.
  
Input parameters are automatically adjusted and results are saved in file 'capon-vertical.max'. By default,
+
The input parameters are automatically adjusted and results are saved in file 'capon-C_135_405-Z.max'. By default,
 
* the minimum velocity is set to 50 m/s,
 
* the minimum velocity is set to 50 m/s,
* time windows are 100-period long,
+
* time windows are 50-period long,
* the minimum frequency is set to have at least 50 time windows,
+
* the cross-spectrum matrix is calculated by averaging 2N blocks (time windows) where N is the the number of sensors,
 +
* the minimum frequency is set to have at least two velocity estimations (taking block average parameters and gaps into account),
 
* the maximum frequency is set to Nyquist frequency,
 
* the maximum frequency is set to Nyquist frequency,
* 2N blocks are used to calculate the cross spectral matrices where N is the number of sensors.
+
* the number of peaks searched on the FK map is limited to N with a minimum relative threshold of 90% for their amplitude (the reference is the highest peak).
  
These values can be changed by providing a parameter file:
+
[[Image:Mirandola-capon-limited.png|thumb|right|500px|Dispersion curve obtained with the vertical component high resolution FK (manual limits).]]
  
  $ geopsy-fk -db Lep_ring01_3C.gpy -group-path vertical -param my.param -o my
+
The results can be viewed with:
  
where ''my.param'' can be for instance:
+
  $ gpviewmax capon-C_135_405-Z.max
  
  PERIOD_COUNT=100
+
The high frequency limit based on Nyquist frequency is far too high for this array. The aliasing limit deduced from Kmax (Wathelet et al., 2008<ref name="Wathelet(2008)">Wathelet, M., D. Jongmans, M. Ohrnberger, and S. Bonnefoy-Claudet (2008). Array performances for ambient vibrations on a shallow structure and consequences over Vs inversion. ''Journal of Seismology'', '''12''', 1-19. , [https://doi.org/10.1007/s10950-007-9067-x 10.1007/s10950-007-9067-x]</ref>) is too restrictive. A curve can certainly be delineated above 4 Hz. Processing parameters can be adjusted in a file:
  MINIMUM_FREQUENCY=1
 
  MAXIMUM_FREQUENCY=30
 
  MIN_V=125
 
  BLOCK_COUNT_FACTOR=1
 
  RELATIVE_THRESHOLD (%)=10
 
  
Parameter ''RELATIVE_THRESHOLD'' selects all FK peaks whose amplitude is higher that 10% of the maximum peak.
+
  $ geopsy-fk -db Mirandola.gpy -group C_135_405-Z -param limits.param -o limited
Option '''-o''' adds a prefix to the automatic output file name which is by default ''process_type-group_name.max''. Alternatively, parameters can be modified directly in the command line:
 
  
  $ geopsy-fk -db Lep_ring01_3C.gpy -group-path vertical -param my.param -o my -set-param "MIN_V=100"
+
where ''limits.param'' can be for instance:
  
In this last case, the parameters are first loaded from ''my.param'' and then ''MIN_V'' is set to 100 instead of 125 m/s. The order of the options '''-param''' and '''-set-param''' matters and they can be used several times.
+
  GRID_SIZE_FACTOR=4
 +
  MAXIMUM_FREQUENCY=10
 +
  MIN_V=150
  
There are many other secondary parameters. To get a list of all possible parameters and their default values, run:
+
[[Image:Mirandola-capon-limited-over100.png|thumb|right|500px|Dispersion curve obtained with the vertical component high resolution FK (manual limits and maximum overlap at 100%).]]
 +
[[Image:Mirandola-capon-limited-over100-stat0.png|thumb|right|500px|Dispersion curve obtained with the vertical component high resolution FK (manual limits, maximum overlap at 100%, no limitation for the number of phase velocity samples).]]
  
  $ geopsy-fk -param-example
+
Option '''-o''' adds a prefix to the automatic output file name which is by default ''process_type-group_name.max''. Parameters can be also modified directly in the command line:
  
[[Image:Lep_dc_vertical_hrfk.png|thumb|right|400px|Dispersion curve for the high resolution FK processed on the vertical component.]]
+
  $ geopsy-fk -db Mirandola.gpy -group C_135_405-Z -param limits.param -o limited -set STATISTIC_MAX_OVERLAP=100 -set STATISTIC_COUNT=0
  
Results can be viewed with:
+
The parameters are first loaded from ''limits.param'', then ''STATISTIC_MAX_OVERLAP'' is set to 100% instead of 0% and ''STATISTIC_COUNT'' is set to 0 instead of 50. The order of the options '''-param''' and '''-set''' matters and they can be used several times. At 100% overlap, successive block sets re-use the same same blocks, sliding by the minimum increment, that is, just one block. The number of phase velocity samples per frequency band is controlled by ''STATISTIC_COUNT''. According to length of blocks and the number of required samples, the block set increment might be larger than the minimum of one. A null ''STATISTIC_COUNT'' removes the limitation of the number of samples (time consuming, particularly at high frequency).
  
  $ gphistogram my-vertical.max
+
There are many other secondary parameters. To get a list of all possible parameters and their default values, run:
  
''.max'' files save the command line arguments, the complete version of ''geopsy-fk'' and its dependencies, the parameters and the results. The history of arguments can be also accessed through
+
  $ geopsy-fk -param-example
 +
 
 +
''.max'' files save the command line arguments, the version of ''geopsy-fk'' and its dependencies, the parameters and the results. The history of arguments can be also accessed through
  
 
   $ geopsy-fk -args
 
   $ geopsy-fk -args
  
Both commands support online help with
+
''geopsy-fk'' and ''gpviewmax'' support help option.
  
 
   $ geopsy-fk -h all
 
   $ geopsy-fk -h all
 +
 +
The next table provides a rough idea of the time spent for various configurations. These estimations were all made on the same machine, an Intel(R) Core(TM) i7-8650U at 1.9 GHz (cache: L1=32 KB, L2=256 KB, L3=8192 KB) with 4 physical cores (+4 with hyper-threading). Two hours of recording with 15 sensors in each case. All runs were performed with release 3.5.2 or 3.5.3 (no difference for FK processing) unless otherwise specified.
 +
 +
{| border="1" cellpadding="5" cellspacing="0"
 +
! Parameters
 +
! Platform
 +
! Threads
 +
! CPU time
 +
! Results
 +
|-
 +
|Default
 +
|Debian 11 (Qt 5)
 +
|4
 +
|58 sec.
 +
|[https://www.geopsy.org/bigfiles/capon-C_135_405-Z.max max]
 +
|-
 +
|Default, direct file access
 +
|Debian 11 (Qt 5)
 +
|4
 +
|1 min. 15 sec.
 +
|[https://www.geopsy.org/bigfiles/capon-importedsignals.max max]
 +
|-
 +
|limits.param, command line
 +
|Debian 11 (Qt 5)
 +
|2
 +
|1 min. 13 sec.
 +
|[https://www.geopsy.org/bigfiles/limited-2jobs-capon-C_135_405-Z.max max]
 +
|-
 +
|limits.param, command line
 +
|Debian 11 (Qt 5)
 +
|4
 +
|46 sec.
 +
|[https://www.geopsy.org/bigfiles/limited-4jobs-capon-C_135_405-Z.max max]
 +
|-
 +
|limits.param, command line
 +
|Debian 11 (Qt 5)
 +
|8
 +
|37 sec.
 +
|[https://www.geopsy.org/bigfiles/limited-8jobs-capon-C_135_405-Z.max max]
 +
|-
 +
|limits.param, graghical interface
 +
|Debian 11 (Qt 5)
 +
|4
 +
|39 sec.
 +
|[[Media:gui-debian11-4jobs-capon-C_135_405-Z.max|max]]
 +
|-
 +
|limits.param, STATISTIC_MAX_OVERLAP=100
 +
|Debian 11 (Qt 5)
 +
|4
 +
|2 min. 11 sec.
 +
|[https://www.geopsy.org/bigfiles/over-capon-C_135_405-Z.max max]
 +
|-
 +
|limits.param, STATISTIC_MAX_OVERLAP=100, STATISTIC_COUNT=0
 +
|Debian 11 (Qt 5)
 +
|4
 +
|21 min. 29 sec.
 +
|[https://www.geopsy.org/bigfiles/allstat-capon-C_135_405-Z.max max]
 +
|-
 +
|limits.param, command line
 +
|VirtualBox Debian 12 (Qt 6)
 +
|4
 +
|34 sec.
 +
|[[Media:debian12-4jobs-capon-C_135_405-Z.max|max]]
 +
|-
 +
|limits.param, command line
 +
|VirtualBox Debian 12 (Qt 6)
 +
|8
 +
|28 sec.
 +
|[[Media:debian12-8jobs-capon-C_135_405-Z.max|max]]
 +
|-
 +
|limits.param, command line
 +
|VirtualBox Window 10 (Qt 6)
 +
|4
 +
|50 sec.
 +
|[[Media:win-4jobs-capon-C_135_405-Z.max|max]]
 +
|-
 +
|limits.param, command line
 +
|VirtualBox Window 10 (Qt 6)
 +
|8
 +
|35 sec.
 +
|[[Media:win-8jobs-capon-C_135_405-Z.max|max]]
 +
|-
 +
|limits.param, graghical interface
 +
|VirtualBox Window 10 (Qt 6)
 +
|4
 +
|50 sec.
 +
|[[Media:gui-win-4jobs-capon-C_135_405-Z.max|max]]
 +
|-
 +
|limits.param, graghical interface, version 3.4.2
 +
|VirtualBox Window 10 (Qt 5)
 +
|8
 +
|3 min. 59 sec.
 +
|[[Media:win-342-capon-C_135_405-Z.max|max]]
 +
|}
  
 
<br style="clear: both"/>
 
<br style="clear: both"/>
Line 117: Line 249:
 
== Single component conventional FK ==
 
== Single component conventional FK ==
  
[[Image:Lep_dc_vertical_convfk.png|thumb|right|400px|Dispersion curve for the conventional FK processed on the vertical component.]]
+
[[Image:Mirandola_conventional.png|thumb|right|500px|Dispersion curve for the conventional FK processed on the vertical component.]]
  
 
To run a conventional FK without block averaging as in ''geopsy'' graphical user interface before 2018:
 
To run a conventional FK without block averaging as in ''geopsy'' graphical user interface before 2018:
  
   $ geopsy-fk -db Lep_ring01_3C.gpy -group-path vertical -param conv.param
+
   $ geopsy-fk -db Mirandola.gpy -group C_135_405-Z -param limits.param -set PROCESS_TYPE=Conventional \
 +
  -set BLOCK_COUNT=1 -set STATISTIC_COUNT=0 -set N_MAXIMA=1
 +
 
 +
{| border="1" cellpadding="5" cellspacing="0"
 +
! Parameters
 +
! Platform
 +
! Threads
 +
! CPU time
 +
! Results
 +
|-
 +
|limits.param, conventional like before 2018
 +
|Debian 11 (Qt 5)
 +
|4
 +
|2 min. 17 sec.
 +
|[https://www.geopsy.org/bigfiles/conventional-C_135_405-Z.max max]
 +
|}
 +
 
 +
<br style="clear: both"/>
 +
 
 +
== All-component ellipticity Rayleigh Direct Steering (ARDS) ==
 +
 
 +
[[Image:Mirandola-ards.png|thumb|right|500px|Dispersion amd ellipticity curves for Rayleigh with ARDS.]]
 +
 
 +
The method is described in Wathelet (2024) <ref name="Wathelet(2024)">Marc Wathelet (2024). Incoherent noise-induced distortions of Rayleigh wave ellipticity measurements obtained with three-component beamforming, Geophysical Journal International, '''236'''(3), 1804-1827, [https://doi.org/10.1093/gji/ggae017 10.1093/gji/ggae017]</ref>.
 +
 
 +
  $ geopsy-fk -db Mirandola.gpy -group C_135_405-3C -param limits.param -set STATISTIC_MAX_OVERLAP=100 -set MINIMUM_FREQUENCY=0.5
 +
 
 +
With only two hours, the default parameters lead to a minimum frequency of 1.27 Hz. Interesting information can be obtained below even if the number of phase velocity samples is lower than STATISTICAL_COUNT.
  
where ''conv.param'' can be for instance:
+
A three-component processing is mush heavier than a single component processing. It sounds reasonable to run it on a HPC infrastructure (high performance computer) rather than a on laptop. A HPC has necessarily a task manager that distribute the cores among the running processes. By default, ''geopsy-fk'' is creating a number of threads equal to half of the available cores (equal to the number of physical cores excluding Hyper-threading). You can adjust manually this number with option '''-j <Nthreads>'''. Alternatively, you can use option '''-task-manager''' which assigns a CPU affinity to each thread. The number of CPU affinities available to your process is controlled by the task manager and indirectly by the number of requested cores for your job. One advantage of assigning a CPU affinity is that each thread is executed by one core during the whole process. On a usual laptop, a thread switches from one core to another several hundreds of times per job. A tiny gain in CPU time can be expected with option '''-task-manager''' on a laptop but only one ''geopsy-fk'' process can by run at time because they will share the same cores rather than being distributed on all cores. Usually it is turned off.
  
   MINIMUM_FREQUENCY=1
+
   $ geopsy-fk -db Mirandola.gpy -group C_135_405-3C -param limits.param -set STATISTIC_MAX_OVERLAP=100 -set MINIMUM_FREQUENCY=0.5 -task-manager
  MAXIMUM_FREQUENCY=30
 
  FREQ_BAND_WIDTH=0.1
 
  PROCESS_TYPE=Conventional
 
  BLOCK_COUNT=1
 
  STATISTIC_COUNT=0
 
  MIN_V (m/s)=100
 
  RELATIVE_THRESHOLD (%)=10
 
  
<br style="clear: both"/>
+
The phase velocity, the ellipticity, the noise parameters, the azimuth and the beam power are view with
 +
 
 +
  $ gpviewmax ards-3C.max
 +
 
 +
{| border="1" cellpadding="5" cellspacing="0"
 +
! Parameters
 +
! Platform
 +
! Threads
 +
! CPU time
 +
! Results
 +
|-
 +
|limits.param, STATISTIC_MAX_OVERLAP=100
 +
|Debian 11 (Qt 5)
 +
|4
 +
|54 min. 29 sec.
 +
|[https://www.geopsy.org/bigfiles/ards-C_135_405-3C.max max]
 +
|-
 +
|limits.param, STATISTIC_MAX_OVERLAP=100, -task-manager
 +
|Debian 11 (Qt 5)
 +
|4
 +
|54 min. 16 sec.
 +
|[https://www.geopsy.org/bigfiles/tm-ards-C_135_405-3C.max max]
 +
|}
  
 
== Rayleigh Three-component BeamForming (RTBF) ==
 
== Rayleigh Three-component BeamForming (RTBF) ==
  
[[Image:Lep_dc_rayleigh_rtbf.png|thumb|right|400px|Dispersion curve for Rayleigh with RTBF.]]
+
[[Image:Mirandola-rtbf.png|thumb|right|500px|Dispersion and ellipticity curves for Rayleigh with RTBF.]]
[[Image:Lep_ell_rayleigh_rtbf.png|thumb|right|400px|Ellipticity curve for Rayleigh with RTBF.]]
+
 
[[Image:Lep_noise_rayleigh_rtbf.png|thumb|right|400px|Estimation of the ratio of incoherent over coherent noise with RTBF.]]
+
The method is described in Wathelet et al. (2018) <ref name="Wathelet et al. (2018)"> Wathelet, M, Guillier, B, Roux, P, Cornou, C. and Ohrnberger, M., Rayleigh wave three-component beamforming: signed ellipticity assessment from high-resolution frequency-wavenumber processing of ambient vibration arrays, ''Geophysical Journal International'', '''215'''(1), 507-523, [https://doi.org/10.1093/gji/ggy286 10.1093/gji/ggy286]</ref>.
[[Image:Lep_dc_love_rtbf.png|thumb|right|400px|Dispersion curve for Love with RTBF.]]
 
  
The method is described in Wathelet et al. (2018) <ref name="Wathelet et al. (2018)"> Wathelet, M, Guillier, B, Roux, P, Cornou, C. and Ohrnberger, M., Rayleigh wave three-component beamforming: signed ellipticity assessment from high-resolution frequency-wavenumber processing of ambient vibration arrays, Geophysical Journal International, 215(1), 507-523.</ref>. A typical parameter file can be (e.g. ''rtbf.param''):
+
  $ geopsy-fk -db Mirandola.gpy -group C_135_405-3C -param limits.param -set PROCESS_TYPE=RTBF -set STATISTIC_MAX_OVERLAP=100 -set MINIMUM_FREQUENCY=0.5
  
  MINIMUM_FREQUENCY=1
+
{| border="1" cellpadding="5" cellspacing="0"
  MAXIMUM_FREQUENCY=30
+
! Parameters
  FREQ_BAND_WIDTH=0.1
+
! Platform
  PROCESS_TYPE=RTBF
+
! Threads
  ROTATE_STEP_COUNT=72
+
! CPU time
  BLOCK_COUNT_FACTOR=4
+
! Results
  MIN_V (m/s)=100
+
|-
  STATISTIC_MAX_OVERLAP(%)=100
+
|limits.param, STATISTIC_MAX_OVERLAP=100
  RELATIVE_THRESHOLD (%)=10
+
|Debian 11 (Qt 5)
 +
|4
 +
|47 min. 45 sec.
 +
|[https://www.geopsy.org/bigfiles/rtbf-C_135_405-3C.max max]
 +
|}
  
Parameter ''STATISTIC_MAX_OVERLAP'' adjusts the overlap between two block sets. At 100%, the block set is shifted by one block. At 0%, the block set is shifted by the number of blocks in the block set to avoid any overlap. Any intermediate value is possible.
+
<br style="clear: both"/>
To run a three-component FK:
 
  
  $ geopsy-fk -db Lep_ring01_3C.gpy -group-path 3C -param rtbf.param
+
== Love beam forming ==
  
Note that the selected group must have the three components available for all sensors.
+
[[Image:Mirandola-capontransverse.png|thumb|right|500px|Dispersion curve for Love with Capon transverse.]]
  
To view the Rayleigh dispersion curve:
+
Love dispersion curve is computed in the same way as in Fäh et al. (2008) <ref name="Fah (2008)">Fäh, D., Stamm, G., & Havenith, H.-B., 2008. Analysis of three-component ambient vibration array measurements, ''Geophysical Journal International'', '''172'''(1), 199–213, [https://doi.org/10.1111/j.1365-246X.2007.03625.x 10.1111/j.1365-246X.2007.03625.x].
 +
</ref>.
  
   $ gphistogram rtbf-3C.max -p R
+
   $ geopsy-fk -db Mirandola.gpy -group C_135_405-3C -param limits.param -set PROCESS_TYPE=CaponTransverse -set STATISTIC_MAX_OVERLAP=100
  
The option '''-p''' defines a pattern to select the result lines of file ''.max'' which contains a ''polarization'' column either Vertical, Rayleigh or Love. To view the Love dispersion curve:
+
A group with only the two horizontal components is sufficient.
  
  $ gphistogram rtbf-3C.max -p L
+
{| border="1" cellpadding="5" cellspacing="0"
 +
! Parameters
 +
! Platform
 +
! Threads
 +
! CPU time
 +
! Results
 +
|-
 +
|limits.param, STATISTIC_MAX_OVERLAP=100
 +
|Debian 11 (Qt 5)
 +
|4
 +
|4 min. 36 sec.
 +
|[https://www.geopsy.org/bigfiles/capontransverse-C_135_405-3C.max max]
 +
|}
  
Love dispersion curve is computed in the same way as in Poggi et al. (2010) <ref name="Poggi et al. (2010)">Poggi, V. and Fäh, D., 2010, Estimating Rayleigh wave particle motion from three-component array analysis of
+
<br style="clear: both"/>
ambient vibrations, Geophysical Journal International, 180(1), 251–267</ref>. To view the Rayleigh ellipticity curve:
 
  
  $ gphistogram rtbf-3C.max -p R -ell-angle
+
== Miscellaneous ==
  
To view the ratio of incoherent over coherent noise:
+
The produced .max files contain an application signature in the header. The signature lists the arguments and the code version (with the Git commit ID). To re-run a case, for instance with new options or a new version, the .max file can be directly executed in a shell (from geopsypack-3.4.2).
  
   $ gphistogram rtbf-3C.max -p R -noise
+
   bash old_results-rtbf-group1.max
  
<br style="clear: both"/>
+
New options can be added.
 +
 
 +
  bash old_results-rtbf-group1.max -o new_results
 +
 
 +
The last option avoids overwriting the original file "old_results-rtbf-group1.max". Even if the original arguments had an option "-o", only the last one is kept.
  
 
== References ==
 
== References ==
 
<references/>
 
<references/>

Latest revision as of 19:41, 17 September 2024

Signal database preparation

This tutorial is based on experimental data recorded in Mirandola (Italy) during InterPACIFIC project. The details of the experimental setup are available in Garofalo et al. (2016) [1].

Download signals AVA in SAC for Mirandola from InterPACIFIC. Decompress the ZIP file. FK processing can access data in two ways:

Database

If you have time, you can build your own database. Alternatively, you can skip this tedious step by downloading Mirandola.gpy. A minor adjustment is still necessary.

Creating a database

To build a database, start by importing the signals in geopsy graphical user interface and by creating a group for each array. There are two options to import signals located in various sub-directories:

 $ geopsy MIRANDOLA_SAC/*/*.sac

or by launching geopsy and importing signals inside the graphical interface with the menu File/Import signals/File pattern. Each sub-directory is a synchronous array. Create a group for each one, eventually with option Default array groups. Set the coordinates of the sensors with the menu Edit/Set receivers. The coordinates are provided in the same archive as the signals (files .geom). Finally save the database to a new .gpy file.

Using an existing database

The section applies only if you intend to use a database created on another machine. A database file contains absolute paths to the signal files. On your computer the original paths certainly do not exist. Paths can be fixed inside geopsy graphical interface but it is also possible to do it through the command line, helpful if you are running on a remote server.

 $ geopsy-fk -db Mirandola.gpy -fix-signal-paths
 ----WARNING--- Opening Database...----
 File '/home/wathelem/Sites/Mirandola/MIRANDOLA_SAC/MIR_C_135_405/MIR_C_135_405_CN01_E.sac' does not exist.
 The path may have changed. Would you like to manually select its new location?
   1. Yes <-- default
   2. No
 ?

Hit enter or answer '1', 'y', 'yes' or just press ENTER.

 Show this message again? [y]/n

Hit enter or answer 'y', 'yes' or just press ENTER.

 ---- Opening Database... ----
   Filter: Signal file (MIR_C_135_405_CN01_E.sac)
   Current directory: /tmp
   Please select file '/home/wathelem/Sites/Mirandola/MIRANDOLA_SAC/MIR_C_135_405/MIR_C_135_405_CN01_E.sac' in its new location.
 File to open:

Provide the absolute or the relative path to the requested file(s). For instance:

 MIRANDOLA_SAC/MIR_C_135_405/MIR_C_135_405_CN01_E.sac

Files having similar paths are automatically translated. At the end the database is saved with your local paths. You can check that all paths are correct by re-starting the same command: you should get no output message. Under Windows, you must use '\' as separator when specifying the new path. It is probably more efficient to do it with the graphical interface under this platform:

 $ geopsy Mirandola.gpy

The list of available groups (name, starting and ending time, duration):

 $ geopsy-fk -db Mirandola.gpy -groups
 geopsy-fk: C_5_15-3C: 20130829100530.000000 20130829110530.000000 (3600 sec.)
 geopsy-fk: C_5_15-Z: 20130829100530.000000 20130829110530.000000 (3600 sec.)
 geopsy-fk: C_15_45-3C: 20130829113830.000000 20130829125330.960000 (4500.96 sec.)
 geopsy-fk: C_15_45-Z: 20130829113830.000000 20130829125330.960000 (4500.96 sec.)
 geopsy-fk: C_45_135-3C: 20130829134300.000000 20130829145600.960000 (4380.96 sec.)
 geopsy-fk: C_45_135-Z: 20130829134300.000000 20130829145600.960000 (4380.96 sec.)
 geopsy-fk: C_26_78-3C: 20130830063630.000000 20130830075631.130000 (4801.13 sec.)
 geopsy-fk: C_26_78-Z: 20130830063630.000000 20130830075631.130000 (4801.13 sec.)
 geopsy-fk: C_135_405-3C: 20130829155630.000000 20130829175500.000000 (7110 sec.)
 geopsy-fk: C_135_405-Z: 20130829155630.000000 20130829175500.000000 (7110 sec.)

Direct file access

This option works only if the coordinates are set in the signal file headers which is not the case for proposed example. The quotes are mandatory to avoid the expansion of file names by the shell.

 $ geopsy-fk "MIRANDOLA_SAC/MIR_C_135_405/*.sac"
 [...]
 geopsy-fk: Relative coordinates with reference at (0.00 0.00).
 geopsy-fk: Null coordinates
 geopsy-fk: Error initializing group 'All signals'

To check the coordinates:

 $ geoosy-fk "MIRANDOLA_SAC/MIR_C_135_405/*.sac" -coord
 0 0 0 CN01
 0 0 0 CN02
 0 0 0 CN03
 [...]

The option -utm can be used to assign coordinates to the signals. The accepted format is relatively rigid: utm_zone x y station_name. The file provided by InterPACIFIC project were edited under Windows, that is, the lines are all ended by CR and LF. Under Linux, a single LF is expected. The following script removes CR and re-orders the columns to fit the expected format.

 cat MIRANDOLA_GEOM/MIR_C_135_405.geom | sed "s/\r$//" |  \
 awk '{if($2!="" && NR>1) {print "---", $2, $3, $1}}' > MIRANDOLA_GEOM/MIR_C_135_405.coord

You can check that the coordinates are properly assigned to the signals.

 $ geopsy-fk "MIRANDOLA_SAC/MIR_C_135_405/*.sac" -utm MIRANDOLA_GEOM/MIR_C_135_405.coord -coord

Single component high resolution FK

Dispersion curve obtained with the vertical component high resolution FK (default parameters).

The simplest way to invoke geopsy-fk is:

 $ geopsy-fk -db Mirandola.gpy -group C_135_405-Z

or

 $ geopsy-fk "MIRANDOLA_SAC/MIR_C_135_405/*.sac" -utm MIRANDOLA_GEOM/MIR_C_135_405.coord -set PROCESS_TYPE=Capon

according to the method chosen in the previous section. In the second option, you must specify the process type to Capon (single component high resolution FK) unless you want a three component process.

The input parameters are automatically adjusted and results are saved in file 'capon-C_135_405-Z.max'. By default,

  • the minimum velocity is set to 50 m/s,
  • time windows are 50-period long,
  • the cross-spectrum matrix is calculated by averaging 2N blocks (time windows) where N is the the number of sensors,
  • the minimum frequency is set to have at least two velocity estimations (taking block average parameters and gaps into account),
  • the maximum frequency is set to Nyquist frequency,
  • the number of peaks searched on the FK map is limited to N with a minimum relative threshold of 90% for their amplitude (the reference is the highest peak).
Dispersion curve obtained with the vertical component high resolution FK (manual limits).

The results can be viewed with:

 $ gpviewmax capon-C_135_405-Z.max

The high frequency limit based on Nyquist frequency is far too high for this array. The aliasing limit deduced from Kmax (Wathelet et al., 2008[2]) is too restrictive. A curve can certainly be delineated above 4 Hz. Processing parameters can be adjusted in a file:

 $ geopsy-fk -db Mirandola.gpy -group C_135_405-Z -param limits.param -o limited

where limits.param can be for instance:

 GRID_SIZE_FACTOR=4
 MAXIMUM_FREQUENCY=10
 MIN_V=150
Dispersion curve obtained with the vertical component high resolution FK (manual limits and maximum overlap at 100%).
Dispersion curve obtained with the vertical component high resolution FK (manual limits, maximum overlap at 100%, no limitation for the number of phase velocity samples).

Option -o adds a prefix to the automatic output file name which is by default process_type-group_name.max. Parameters can be also modified directly in the command line:

 $ geopsy-fk -db Mirandola.gpy -group C_135_405-Z -param limits.param -o limited -set STATISTIC_MAX_OVERLAP=100 -set STATISTIC_COUNT=0

The parameters are first loaded from limits.param, then STATISTIC_MAX_OVERLAP is set to 100% instead of 0% and STATISTIC_COUNT is set to 0 instead of 50. The order of the options -param and -set matters and they can be used several times. At 100% overlap, successive block sets re-use the same same blocks, sliding by the minimum increment, that is, just one block. The number of phase velocity samples per frequency band is controlled by STATISTIC_COUNT. According to length of blocks and the number of required samples, the block set increment might be larger than the minimum of one. A null STATISTIC_COUNT removes the limitation of the number of samples (time consuming, particularly at high frequency).

There are many other secondary parameters. To get a list of all possible parameters and their default values, run:

 $ geopsy-fk -param-example

.max files save the command line arguments, the version of geopsy-fk and its dependencies, the parameters and the results. The history of arguments can be also accessed through

 $ geopsy-fk -args

geopsy-fk and gpviewmax support help option.

 $ geopsy-fk -h all

The next table provides a rough idea of the time spent for various configurations. These estimations were all made on the same machine, an Intel(R) Core(TM) i7-8650U at 1.9 GHz (cache: L1=32 KB, L2=256 KB, L3=8192 KB) with 4 physical cores (+4 with hyper-threading). Two hours of recording with 15 sensors in each case. All runs were performed with release 3.5.2 or 3.5.3 (no difference for FK processing) unless otherwise specified.

Parameters Platform Threads CPU time Results
Default Debian 11 (Qt 5) 4 58 sec. max
Default, direct file access Debian 11 (Qt 5) 4 1 min. 15 sec. max
limits.param, command line Debian 11 (Qt 5) 2 1 min. 13 sec. max
limits.param, command line Debian 11 (Qt 5) 4 46 sec. max
limits.param, command line Debian 11 (Qt 5) 8 37 sec. max
limits.param, graghical interface Debian 11 (Qt 5) 4 39 sec. max
limits.param, STATISTIC_MAX_OVERLAP=100 Debian 11 (Qt 5) 4 2 min. 11 sec. max
limits.param, STATISTIC_MAX_OVERLAP=100, STATISTIC_COUNT=0 Debian 11 (Qt 5) 4 21 min. 29 sec. max
limits.param, command line VirtualBox Debian 12 (Qt 6) 4 34 sec. max
limits.param, command line VirtualBox Debian 12 (Qt 6) 8 28 sec. max
limits.param, command line VirtualBox Window 10 (Qt 6) 4 50 sec. max
limits.param, command line VirtualBox Window 10 (Qt 6) 8 35 sec. max
limits.param, graghical interface VirtualBox Window 10 (Qt 6) 4 50 sec. max
limits.param, graghical interface, version 3.4.2 VirtualBox Window 10 (Qt 5) 8 3 min. 59 sec. max


Single component conventional FK

Dispersion curve for the conventional FK processed on the vertical component.

To run a conventional FK without block averaging as in geopsy graphical user interface before 2018:

 $ geopsy-fk -db Mirandola.gpy -group C_135_405-Z -param limits.param -set PROCESS_TYPE=Conventional \
 -set BLOCK_COUNT=1 -set STATISTIC_COUNT=0 -set N_MAXIMA=1
Parameters Platform Threads CPU time Results
limits.param, conventional like before 2018 Debian 11 (Qt 5) 4 2 min. 17 sec. max


All-component ellipticity Rayleigh Direct Steering (ARDS)

Dispersion amd ellipticity curves for Rayleigh with ARDS.

The method is described in Wathelet (2024) [3].

 $ geopsy-fk -db Mirandola.gpy -group C_135_405-3C -param limits.param -set STATISTIC_MAX_OVERLAP=100 -set MINIMUM_FREQUENCY=0.5

With only two hours, the default parameters lead to a minimum frequency of 1.27 Hz. Interesting information can be obtained below even if the number of phase velocity samples is lower than STATISTICAL_COUNT.

A three-component processing is mush heavier than a single component processing. It sounds reasonable to run it on a HPC infrastructure (high performance computer) rather than a on laptop. A HPC has necessarily a task manager that distribute the cores among the running processes. By default, geopsy-fk is creating a number of threads equal to half of the available cores (equal to the number of physical cores excluding Hyper-threading). You can adjust manually this number with option -j <Nthreads>. Alternatively, you can use option -task-manager which assigns a CPU affinity to each thread. The number of CPU affinities available to your process is controlled by the task manager and indirectly by the number of requested cores for your job. One advantage of assigning a CPU affinity is that each thread is executed by one core during the whole process. On a usual laptop, a thread switches from one core to another several hundreds of times per job. A tiny gain in CPU time can be expected with option -task-manager on a laptop but only one geopsy-fk process can by run at time because they will share the same cores rather than being distributed on all cores. Usually it is turned off.

 $ geopsy-fk -db Mirandola.gpy -group C_135_405-3C -param limits.param -set STATISTIC_MAX_OVERLAP=100 -set MINIMUM_FREQUENCY=0.5 -task-manager

The phase velocity, the ellipticity, the noise parameters, the azimuth and the beam power are view with

 $ gpviewmax ards-3C.max
Parameters Platform Threads CPU time Results
limits.param, STATISTIC_MAX_OVERLAP=100 Debian 11 (Qt 5) 4 54 min. 29 sec. max
limits.param, STATISTIC_MAX_OVERLAP=100, -task-manager Debian 11 (Qt 5) 4 54 min. 16 sec. max

Rayleigh Three-component BeamForming (RTBF)

Dispersion and ellipticity curves for Rayleigh with RTBF.

The method is described in Wathelet et al. (2018) [4].

 $ geopsy-fk -db Mirandola.gpy -group C_135_405-3C -param limits.param -set PROCESS_TYPE=RTBF -set STATISTIC_MAX_OVERLAP=100 -set MINIMUM_FREQUENCY=0.5
Parameters Platform Threads CPU time Results
limits.param, STATISTIC_MAX_OVERLAP=100 Debian 11 (Qt 5) 4 47 min. 45 sec. max


Love beam forming

Dispersion curve for Love with Capon transverse.

Love dispersion curve is computed in the same way as in Fäh et al. (2008) [5].

 $ geopsy-fk -db Mirandola.gpy -group C_135_405-3C -param limits.param -set PROCESS_TYPE=CaponTransverse -set STATISTIC_MAX_OVERLAP=100

A group with only the two horizontal components is sufficient.

Parameters Platform Threads CPU time Results
limits.param, STATISTIC_MAX_OVERLAP=100 Debian 11 (Qt 5) 4 4 min. 36 sec. max


Miscellaneous

The produced .max files contain an application signature in the header. The signature lists the arguments and the code version (with the Git commit ID). To re-run a case, for instance with new options or a new version, the .max file can be directly executed in a shell (from geopsypack-3.4.2).

 bash old_results-rtbf-group1.max

New options can be added.

  bash old_results-rtbf-group1.max -o new_results

The last option avoids overwriting the original file "old_results-rtbf-group1.max". Even if the original arguments had an option "-o", only the last one is kept.

References

  1. Garofalo, F., Foti, S., Hollender, F., Bard, P., Cornou, C., Cox, B. R., Ohrnberger, M., Sicilia, D., Asten, M., Di Giulio, G., Forbriger, T., Guillier, B., Hayashi, K., Martin, A., Matsushima, S., Mercerat, D., Poggi, V., Yamanaka, H. (2016). InterPACIFIC project: comparison of invasive and non-invasive methods for seismic site characterization. Part I: intra-comparison of surface wave methods, Soil Dynamics and Earthquake Engineering, 82, 222-240, 10.1016/j.soildyn.2015.12.010
  2. Wathelet, M., D. Jongmans, M. Ohrnberger, and S. Bonnefoy-Claudet (2008). Array performances for ambient vibrations on a shallow structure and consequences over Vs inversion. Journal of Seismology, 12, 1-19. , 10.1007/s10950-007-9067-x
  3. Marc Wathelet (2024). Incoherent noise-induced distortions of Rayleigh wave ellipticity measurements obtained with three-component beamforming, Geophysical Journal International, 236(3), 1804-1827, 10.1093/gji/ggae017
  4. Wathelet, M, Guillier, B, Roux, P, Cornou, C. and Ohrnberger, M., Rayleigh wave three-component beamforming: signed ellipticity assessment from high-resolution frequency-wavenumber processing of ambient vibration arrays, Geophysical Journal International, 215(1), 507-523, 10.1093/gji/ggy286
  5. Fäh, D., Stamm, G., & Havenith, H.-B., 2008. Analysis of three-component ambient vibration array measurements, Geophysical Journal International, 172(1), 199–213, 10.1111/j.1365-246X.2007.03625.x.