From GeopsyWiki
Jump to: navigation, search

F-K Toolbox (conventional f-k) - overview

Signal viewer for one group of simultaneously recorded waveforms ready for array processing

If you have followed all steps for the preparation of a database, loading coordinates and grouping, you should have now a set of simultaneously recorded waveforms building a group and containing coordinates for each waveform.

Alternatively, you can download here a ready-to-go database file and the corresponding waveform files at this location. Note that you have to re-locate your waveform files as described in detail in this page. This database contains two predefined groups: vertical and All 3C. Please test, whether you can view the signals, use drag and drop functionality or check the coordinate settings in the table view or map view. If you display the predefined group vertical in a graphic viewer, you should obtain a picture similar to the one displayed on the right.

Now you can open the conventional frequency wavenumber toolbox by pushing the plugin icon FKToolboxPluginIcon.png.

The f-k toolbox should open and is now attached to the 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. Note, for reasons understandable from section testing parameter settings it will be always a graphic signal viewer that is openend.

Conventional f-k toolbox

Let's now have a close look to the individual parameters to be set in the f-k toolbox. There are four main groups of parameters to be considered.

  • Pre-processing parameters for excluding data from processing.
  • Processing parameters for selecting time windows on the remaining data.
  • Processing parameters for selecting narrow frequency bands for processing
  • Parameters related to the output of results from f-k processing.

These parameters are laid out in different parts of the f-k toolbox widget. Note, that the toolbox has three main tabs: Time, Processing and Status.

The pre-processing parameters for excluding data from processing follow the same argumentation as the data selection for other toolboxes like H/V processing and can be found in the Time tab.

E.g. one can apply the so-called anti-trigger 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 anti-trigger. Another, second, approach for excluding data from processing is provided via the bad sample approach. For details in setting parameters with the bad sample approach check this page.

Parameter settings

Parameters for time window settings

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.

For details on possible options for setting the time limits, please have a look here

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 FKWindowLength.png.

  • 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 reasonable and only implemented setting for the Length parameter in the context of narrowband dispersion curve processing is Freq. dep.. For the conventional frequency wavenumber values range typically between 20.00 to 50.00 meaning that 20 to 50 cycles of the central signal 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. In upcoming versions of geopsy software, this combo box entries will be adjusted to the reasonable options.

Parameters for frequency band selection

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 FKFreqSampling.png.

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 FKGridding.png.

The frequency limits for each frequency band are computed as 
[(1-bw)\cdot f_c,(1+bw)\cdot f_c]
where f_c is the central frequency and bw 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:

Frequency band creation using Step: Linear, Bandwidth: 0.10, From: 1 Hz, to: 10 Hz
Freq. band index 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 \vec{k}=(k_x,k_y)^T and computing the semblance coefficient Semb(\omega,\vec{k}) (coherence measure, see Neiddell and Taner (1971) [1] or Douze and Laster (1979) [2]) and/or beam power BP(\omega,\vec{k}) of the shifted stacked output of all array stations. Beampower and Semblance are defined as follows:

  • Beampower:

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

  • Semblance:

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}
{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}

In the above formulas X_i(\omega) is the observed record at station i at frequency \omega. The station coordinates in the array plane (usually earth's surface) are given by \vec{r}_i. Therefore, the time delays of a harmonic plane wave propagating along the horizontal plane in direction of \vec{k} and apparent propagation velocity c related to the absolute length of the wavenumber vector \|\vec{k}\| = 2\pi/\lambda = 2\pi f/c = \omega/c can be computed as \tau = \vec{k}\vec{r}_i. Thus, the beampower/semblance measure the power/coherence of a plane wave propagating with wavenumber vector \vec{k}. By testing many wavenumber vectors in the wavenumber plane one tries to find those wavenumber vectors 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 (k_x,k_y)^T). In order to optimize the grid search approach, the f-k toolbox relies on an iterative refined sampling strategy. In the first iteration potential maxima are identified within the search limits of the wavenumber plane using a coarse grid step in both axis directions (k_x,k_y). Note, that coarse is relative and must be fine enough in order to match the goal: identifying potential maxima in the wavenumber plane. Therefore, the grid step must be as small enough as to sample the peaky structure of the array reponse sufficiently fine for not missing any small maximum in the optimization procedure. The array response is closely related to the 2-dimensional array geometry (please check also the tools gpfksimulator and warangps for this issue). In particular the width of the main lobe peak - we refer to it as k_{min} - in the array response (which is the one we are looking for) is approximately related to the overall aperture of the array D_{max} by k_{min} \approx 2\pi / D_{max}.

In the current version of the f-k toolbox, the grid_step for sampling is automatically calculated as grid\_step = k_{min}/4. Thus, any potential peak to be identified in the wavenumber domain will be sampled in either direction at least with four samples.

In a similar way, the maximum search radius in the wavenumber space (we call it grid_size) is related on one hand to the minimal wave speed that is expected to be observable in the wave field. The second limit to be considered is related to the existence of strong side lobes or grating lobes in the array reponse function. However, this limit is much more ambiguous than the grid_step. Taking the viewpoint of scanning the whole wavenumber space up to physical limits of wave propagation velocities will lead often to an unnecessary computation time for lower frequency bands. On the other hand, restricting the grid_size too strongly on the apparent side lobes in the array reponse function will not allow for high frequencies to obtain any result for slowly propagating waves. For an in-depth discussion of this issue look here.

Gridding parameters for the wavenumber plane

The setting for the parameters regarding the gridding of the wavenumber plane are summarized in the FK gridding layout group of the processing tab in the f-k toolbox (see picture to the left). The values for grid_step and grid_size are computed automatically following the comments provided about the minimum physical propagation velocity is set per default to 100 m/s. Depending on the knowledge about your site conditions, you may or may not modify this setting. Be aware that soft soil sites may find Rayleigh wave propagation velocities for the fundamental mode fundamental mode well below 100 m/s!

Parameter settings for output / saving results

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 FKLoadParams.png.

Note: the option Maximum number allows you to store more than the highest maximum found in the optimized grid search in the wavenumber plane. In order to identify dispersion curve branches at higher frequencies where a lot of aliasing may occur, it may prove beneficial to use this option (e.g. by replacing the defualt number 1 to a value of 2 to 4).

There are two additional output filtering options that can be used to store not every result into the output files. The options are named Absolute th. and Relative th. and are explained in detail in max2curve. The default values are 0 in both cases meaning that no filtering on the output results will occur. As one is able to filter the result output at a later stage of the processing chain using max2curve we certainly recommend not to change these values. As it is a post-processing filtering there will be also no gain in processing speed.

Testing parameter settings and starting f-k processing

Test button for checking parameter settings and controlling data quality interactively

Once you finished all parameter settings as described above, you may check the impact of the settings on the f-k processing by using the Test button in the Time tab of the f-k toolbox.

interactive frequency wavenumber window browser

Pushing the Test button will open a frequency wavenumber window browser which allows you to explore the f-k analysis for individual time windows and different frequency bands.

You can move the Time window slider to scroll in time over the selected data and selecting the individually processed data windows. The results of the wavenumber analysis will be displayed in the browser window as a colored map. The color scale is following the physical rainbow light spectrum, associating violet with high and red with low power values. While sliding through time, you will notice that the wavenumber results will change depending on the direction of wave propagation and number of identifiable plane wave sources in the wave field.

Note: please be aware of caveats for the interpretation of multiple plane waves propagating at the same time with distinct propagation properties (check the page on multiple plane wave arrivals and explore this issue using gpfksimulator).

An interesting feature of this wavenumber browser window is the capability to control the consistency of apparent wave velocity / horizontal slowness estimation. The black circle, which is displayed on top of the wavenumber plane corresponds to a single propagation velocity / slowness without considering the propagation azimuth. For a wave field completely composed of a mixture of a single mode surface wave (e.g. Rayleigh for vertical component analysis), we would therefore expect the peak in the wavenumber domain always to be located at the same distance from the center of the wavenumber map. Using the black circles helps to locate the relative distance of dominant peak(s) in the wavenumber domain with respect to the center of the figure. For a single frequency, we expect only the direction to change, but finding the peak maximum always below the circle line.

Status tab of the fk toolbox while running an analysis

Now as you have tested your parameter configuration and had a preliminary look on the data, you can start the processing by simply pushing the Start button in the lower right corner of the f-k toolbox (visible in all tabs). You can check the progress of the computation in the status tab of the toolbox. Note that geopsy and its plugins are multi-threaded applications and will use all available cpu-cores of your hardware. The progress for each computation thread will be displayed separately.

Contents of output file

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 entries are obsolete but are kept for compatibility with former file formats used in other codes during the SESAME project.

# File generated by Geopsy, FK processing
# The process log is saved in file lep_ring01.log
# Header 
# NPHI 72
# 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 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 ( 10\cdot \log(BP(\omega,\vec{k}))

The Log file contents as displayed below contain three different sections, Init. Parameter and Process. The init section reports which signals of which stations are included for processing and the relative coordinates in units of m are reported for each station. The parameter section contains a summary of all settings specified in the f-k toolbox setting and which are used for the processing. The process section finally contains the summary of all windows processed for each frequency bands. Examples for the three sections are given below:

### 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: when loading the log file via the Load Params button, those parameters will be loaded into the f-k toolbox.

### Parameters ###
WINDOW LENGTH TYPE (at least/exactly/freq. dep.) = freq. dep.
USED RAW COMPONENTS = y, n, y, y, y, y, y, y, y, y
RAW STA (s) = 1
RAW LTA (s) = 30
SAMPLING TYPE FREQUENCY (0=log, 1=linear)= 0
FROM TIME TEXT = 8h21m0s
TO TIME TEXT = 9h29m0s
MIN K = 0.015
MAX K = 0.5
MIN V = 100
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

For the interpretation of the results of the f-k processing one can use the graphical tool max2curve or a combination of gphistogram and max2curve.

You may call max2curve on the command line directly with the output file (.max file) produced from your f-k processing.

$> max2curve your_output_file.max
Startup of max2curve for selecting histogram parameters

Certainly you may also open max2curve directly. Then you will need to locate you .max file in the directory tree and select it with the automatically opened file browser.

The utility program max2curve then reads the contents of the f-k processing output file and computes a probability density function (pdf) for slowness values for each individual frequency band. In order to be able to produce the pdf from re-normalization of a binned histogram, max2curve needs to know about the slowness limits and the number of bins to use. As you have already looked at the picture on the left you will say: stop - you probably mean an apparent velocity histogram or pdf of apparent velocity values for a given frequency.

Smart thought, because the entries for the limits (Minimum velocity and Maximum velocity) are to be specified in (apparent) velocity. And additionally one has to specify the Number of velocity classes. Despite these confusing entries, the computation indeed takes place in the slowness domain with a linear division of the slowness axis into the number of bins chosen.

As you are probably now fully confused, let me explain the thoughts behind it: the computation takes place in slowness domain, because our estimates are obtained from measuring time delays. Thus, the error of this measurement goes linear with the time delay estimates and therefore is linear in slowness. So, in terms of error propagation the computation makes clearly sense in slowness domain and not(!) in the inversely related velocity domain.

So why, you may ask, give the limiting values in apparent velocity anyway? Well, this is a concession to the (still) large part of the community who this used to think in wave velocities and not the inverse property.

As already mentioned above in section parameter settings for output, it is possible to filter the result output (.max file(s)) according to two extra criteria. Using the two entries under Sample selection, i.e. specify percentage values for Semblance and Beam power, it is possible to filter individual time window results falling below thresholds for the semblance coefficients and beam power values. Although semblance is a normalized quantity, we choose to develop a relative threshold scheme for filtering results for both semblance and beam power values.

The percentage values should be understood as follows: for all entries in the .max output file, the minimum and maximum values for the semblance coefficient and the beam power values are used as lower (0%) and upper (100%) limits. Any given percentage value will compute the threshold as absolute value in between these two limits by linear interpolation. Note that the conditions are combined in the sense of a logical and. Only time window that surpass the threshold for both values (semblance and beam power) will be kept and passed to the histogram computation. Any other analysis window result will be skipped and is not taken into account for the computation of the frequency dependent histogram.

max2curve main window after startup

After having set all necessary values (default values may be a good starting point :-)), max2curve will start up with its main window consisting of two main parts. On the left, there is a curve browser element. On the right, the slowness histograms (pdfs) are displayed as a color image for a log frequency scaling. Per default behaviour max2curve computes the sample mean and sample standard deviation for each frequency histogram. Those values are stored in the curve browser as automatically determined mean curve with error estimates.

Note: Due to the variability in the estimation process, various disturbances in the wavefield and violations of the basic array processing assumptions (e.g. plane wave arrivals) and the well-known limitations of the array geometry in terms of aliasing and resolving capabilities, the automatically determined dispersion curves look often very messy. In order not to get disturbed by this curve plotted over the histogram display, you can hide the curve by toggling the visible button an the bottom of the curve browser.

Grid statistics window of max2curve

The Grid statistics window allows the user of max2curve to visualize the histograms for each frequency. To scan over different frequency bands, one can use the slider bar to scan from the lowest to the highest frequency band. For some selected center frequency, one obtains the re-normalized histogram ( = pdf) of all slowness estimates in the upper panel. The sample mean \mu(f_c) and sample standard deviation \sigma(f_c) calculated from the histogram are used to display the corresponding normal distribution N(\mu(f_c),\sigma(f_c))) as overlaying red curve.

The lower panel shows the same information, but displayed as cumulative histogram or probability function instead of histogram and pdf. The display has been introduced because it is often easier to recognize the proximity of a distribution to the normal distribution in the cumulative view.

Besides the visualization, the grid statistics window can be further used to modify/edit/clean the histograms on grounds of expert knowledge about the way aliasing phenomena are wrapped into the histogram picture. Details about the editing capabilities and final dispersion curve determination may be found in this page.

Note on three-component f-k processing

The processing of all three components of an array recording is easy to achieve technically. Instead of using single component data in one group, just define a group containing all three components for each station. In the example database provided above in section overview you would like to use the group named All 3C. Opening then the f-k toolbox for this group will automatically enable the three component f-k processing. For each propagation direction, the horizontal component signals are rotated into radial and transverse directions and the processing will be performed on vertical, radial and transverse components. For the output, three distinct output files will be created, by adding vertical, radial and transverse to the base filename and adding the extensions .max and .log as before.

Please note that the interpretation of the results for the horizontal components may present more difficulties than expected.

The wavenumber window browser allows to view the wavenumber planes individually for vertically, radial and transverse directions.


  1. Neiddell, N. and Turhan Taner, M.: Semblance and other coherency measures for multichannel data, Geophysics, 36, 482–497, 1971.
  2. Douze, E. J. and Laster, S. J.: Statistics of semblance, Geophysics, 44, 1999–2003, doi:10.1190/1.1440953, short Note, 1979.