/ CEA/DASE
PMCC Documentation
Main Features and tuning parameters

1.pmcc: main features......

1.1Introduction......

1.2Detection by correlation......

1.3Time correction for a non planar array......

1.4Threshold of detection: the consistency

1.5The progressivity

1.6Checking data quality

1.7Time frequency analysis and post-processing......

2description and Tuning of the Main Parameters......

2.1Array parameters......

2.2Detection parameters......

2.3Frequency parameters......

2.4Family parameters......

3REFERENCES......

4APPENDIX......

1.pmcc: main features

1.1Introduction

In contrast to a set of isolated sensors, a dense array, whose aperture is of the order of the wavelengths of the signals of interest, allows similarity measurements of the recordings to avoid uncertainties encountered with individual arrival-time picking.

Similarly to seismology, most of the infrasonic waves can be represented at a local scale by a set of planar waves using the well-known relation

where:

is the wave vector associated to frequency f and phase velocity c

is the angular frequency.

The frequency content of a recorded wave can easily be determined using a single station. At the opposite, a set of sensors is needed to calculate the propagation parameter. When the aperture of this set of sensors is equivalent to some wavelengths of signal, this set is named an array. At the opposite, when the aperture is much larger than the wavelength, it is named a network.

In the case of a network, the signal is often very different from one sensor to another and the measure of the propagation parameters is derived from the set of arrival-times by inversion, as described by Husebye (1969) in seismology. On the opposite, in the case of an array, we use the similarity of the signals to compute arrival time differences using classical techniques of signal processing theory. This set of arrival time differences is used to compute the propagation parameters with a Husebye’s derived method.

The most classical method for estimating these wave parameters in the case of an array is a systematic search in a specific domain of wave vector using the signals recorded on the sensors. For example the disc defined in the wavenumber plane by the relation corresponds to all the waves with a frequency f, with any azimuth and with a velocity. For each discrete wave vector of this regularly discretized domain, the time delay at each sensor is calculated and the delayed signals are summed. When the signals are mainly composed of random background noise, the energy variation of the sum is small over the entire wave vector field. In contrast, if the signals are associated with a specified vector, the energy obtained for will be much larger than for the other vectors. A lot of methods have been proposed by different authors to find the wave vector which produces the maximum energy [Capon, 1969]. This is not a trivial problem because data are discrete in the space domain (i.e. only few sensors are used). This implies that for each frequency, false results can be obtained due to correlated signals over one or more periods (ambiguity effect). To study these effects, Capon suggested to compute the array beam forming function. It represents the array response for a wave which arrives vertically under it, i.e. with a horizontal planar wave reaching all the sensors at the same time. The phase velocity of the wave measured is then infinite, which leads to a maximum amplitude for.

The main assumption linked to this type of methods is the search of a signal modeled by a planar wave recorded on all the sensors of the array. In practice, some of these assumptions are confirmed to various extents as diverse sources cause interference with these signals. That is why a more flexible method, less constraining with respect to the model, is proposed. It is based on conventional signal processing techniques to detect a stable signal on two or more records, partly by relaxing the planar wave model rigidity. The PMCC method (Progressive Multi-Channel Correlation), originally designed for seismic arrays, proved also to be efficient for analyzing low-amplitude infrasonic coherent waves within non-coherent noise.

1.2Detection by correlation

A temporal signal can be represented by its Fourier transform where represents the spectral amplitude and is the phase. The background noise is characterized by a rapid variation of both and from one sensor to another, even if they are closer than one wavelength of signal. On the opposite, in case of signal propagating between the sensors, the following relations are available:

These relations indicate that no deformation exists between the two signals, and that the only difference is a delay depending on the relative positions of the sensors (i.e.: in the case of a planar wave).

Based on these two observations, a signal-processing tool can be used to detect a signal present on the records and. The correlation function is used to measure the time delay between the two records. In case of a wave propagating without distortion (assuming a planar wave), this delay is the same for all frequencies of the signals:

This measurement is made in the time domain with values ranging from -1 to 1. Taking into account all frequencies, it measures in a given time window the similarity of the signals when shifted in time. The maximum of the correlation function gives the time delay between the signals. This method enables a decision to be made on whether there is a signal in a set of simultaneous records, independently of any information on previous records.

1.3Time correction for a non planar array

In some cases, for a non planar array, the travel time differences due to elevation differences between sensors become not negligible. For such arrays, the elevation of each array element should be taken into account. The calculation of the time delays between two sensors is to be improved. It can be shown that the time delays then also depends on the local sound speed and the elevation angle. In a first order, this correction is applied:

where:

 is the backazimuth (often called azimuth), angle of wave front approach, measured clockwise between the north and the direction towards the source,

c is the local sound speed

d is the horizontal distance between sensor j and the barycenter

j is the time-delay between sensor j of coordinates (x,y,z) and the center of the array (0,0,0)

i is the incidence angle between the direction of propagation of the wave front and the vertical

, is the horizontal trace velocity across the array.

1.4Threshold of detection: the consistency

To avoid ambiguity problems when correlating the records from sensors too far apart, the analysis is initialized on the smallest groups of three sensors. The correlation function is used to calculate the propagation time tij of the wave between sensors i and j. For each subnetwork (i,j,k), the closure relation tij +tjj+tki = 0 should be obtained. In the presence of background noise the phase is unstable. Therefore, the delays measured in this case are the result of random phase combinations. These delays, independent of the amplitude of each elementary wave, become random, and the closure relation given above is no longer valid.

The consistency of the set of delays obtained using all the sensors of is then defined as the mean quadratic residual of the closure relations:

If this consistency is below a given threshold, a detection is observed on .

1.5The progressivity

To minimize errors in the calculation of the wave parameters, distant sensors are progressively added using a criterion based on a comparison between their distance to the subnetwork and the computed wavelength. This progressive use of distant sensors has two main effects:the removal of false detections which could be due to correlated noise at the scale of the starting subarrays, anda better estimation of the wave parameters by increasing the array aperture.

After being initialized with a small subnetwork of three sensors, in order to avoid ambiguity problems inherent in the correlation of signals from distant sensors, the wave parameters calculated on the initial subnetwork are used when adding other sensors. For that, a propagation of a planar wavefront is assumed. The new measured time delay is given by the maximum of the correlation function which is the closest to the one that has been estimated. Each elementary detection is therefore defined by several parameters such as the consistency value, the number of sensors participating to the detection, the frequency, the horizontal trace velocity and the backazimuth.

Such a detector is independent of the signal amplitude and uses only the intrinsic information of the recordings. As long as the closure relation is valid, the use of sensors increasingly further apart gives more precise wave parameters since the aperture of the network increases with each new sensor. The final solution is given by the largest subnetwork.

1.6Checking data quality

In the PMCC configuration file (see in Appendix), a list of starting subnetworks of three sensors is defined. As explained previously, the detection is initiated through each of these subnetworks. Then, as long as there are available sensors and the consistency value remains lower than its threshold, new sensors are aggregated. To avoid wrong results due to the lack of data in the recordings, an automatic procedure checks the data quality. If the initial subnetworks contain sensors with consecutive zeros in the recordings, this procedure looks for other set of three sensors belonging to the array. Among all possible combinations calculated from the remaining sensors, the best subnetworks are selected. The principle is to sort them according to symmetry and size criteria. Equilateral triangle of small aperture is the best configuration. The maximum number of new eligible subnetworks corresponds to the number of subnetworks defined in the configuration file.

1.7Time frequency analysis and post-processing

The processing is performed consecutively in several frequency bands fi, and in adjacent time windows tj covering the whole period of analysis. To avoid unrealistic wave detection, a further condition is introduced. A set of several elementary detections in the time-frequency domain is considered to represent one detected wave (corresponding for example to different frequency bands or adjacent windows). Conversely, several waves with different parameters may coexist in the same time window but in different frequency bands. Each wave must be identified separately. To do this, a nearest-neighbor search of elementary detections in the (ti,fi,VTi,i) domain is used [Cansi, 1995; Cansi, 1997]. The final detection is thus an aggregate of close-enough points in this domain. Different distances can be used to connect close-enough points, but the more usual is a weighted Euclidian distance defined by:

The ’s are used as weighting factors to allow the comparison of quantities of different meanings. These weights can be tuned independently. As opposed to the other weights, is the only dimensionless parameter. The velocity weight is expressed in percent. Figure 1 presents an example of PMCC results before and after this processing.

Figure 1 - Results of PMCC calculation on signals generated by a gaz explosion in Belgium recorded at the Flers experimental infrasound station set up in Normandy (France). The PMCC results (horizontal trace velocity and azimuth) are presented in time / frequency diagrams. Values are given according to the color scales. The results are presented from 0.1 to 4 Hz in 10 equally spaced frequency bands. Azimuths are given clockwise from North. Top: elementary detection. Bottom: final results after the post-processing,

2description and Tuning of the Main Parameters

Depending on type of signals to be processed (eg waveforms and frequency content) and the background noise level, the PMCC parameters have to be set in order to ensure an optimum detection providing a compromise between a low detection threshold and low numbers of false detections. Too large threshold values will allow a low level of detection, but a large number of false detections are expected. On the opposite, too small values are not appropriate to detect coherent waves with poor signal-to-noise ratio, but a very few number of false detections will be obtained.

PMCC requires input of two configuration files: pmcc.ini and filters.ini. These contain, respectively, the computation parameters and the filter coefficients. Examples of these two files as well as the format of the result files are provided in the Appendix. The main PMCC parameters are related to the array, the detection, the filtering and the post processing. They are described below along with some guidelines underlying the choice of suitable values.

2.1Array parameters

  1. [X,Y] are the sensor coordinates (relative to the barycenter, in km). For a relevant calculation of the wave parameters, they should be accurate enough regarding the sampling rate, the array aperture and the wavelength of the signals.
  2. NbSensors: number of sensors of the arrray.
  3. NbSubNet is number of subnetworks to be processed. In order to avoid ambiguity problems inherent in the correlation of signals from distant sensors, the calculation is generally initiated with small subnetworks of three sensors. Then, for each subnetwork, distant sensors are progressively added using a criterion based on a comparison between their distance to the subnetwork and the computed wavelength. This progressive use of distant sensors has two main effects:the removal of false detections which could be due to correlated noise at the scale of the starting subarrays, anda better estimation of the wave parameters by increasing the array aperture. Since all sensors of the array will be added as long as the consistency remains a threshold value, it is not necessary to define all available subnetworks. Furthermore, the computation time grows linearly with NbSubNet. For an 8-element array, a maximum reasonnable number of NbSubNet is 5. If the aperture of the full-array is significantly smaller than the wavelength, in order to make the process less sensitive to highly correlated undesirable signals, the largest subnetworks can be defined.
  4. GeoSubNetdescribes each selected subnetwork, with numbers referring to the sensors. Their aperture should be in the same order of wavelength of the signals and their geometry should be as symetric as possible. The equilateral triangle is the optimum configuration.
  5. DelayCorrection (s) is the time correction to be applied to each channel (null is the default value).
  6. Detection parameters
  1. WindowLength (s) is the length of time window to compute the correlation. It can not be lower than the maximum propagation time of the wave between two sensors. This value is generally increased to ensure a smoothing of the time base, which significantly reduces the number of false detections and improves statistics on the measurements of the wave parameters. A reduction in the number of families (final detection of PMCC aggregating close-enough pixels in the time / frequency / speed / azimuth domain) is then expected, providing a synthesis of information. In the case of long duration phenomena such as microbaroms or Mountain-Associated Waves, WindowLength can contain several periods (up to 10). In case of transient signals, WindowLength can also be greater than the maximum propagation time between sensors. The disadvantage of a broad time window is the anticipation of the detection and the increase of its duration. In order to take into account this effect, a time picking procedure has to be used for an accurate arrival time measurement.
  2. TimeStep (s)is the time shift between two consecutive windows. It cannot exceed WindowLength. It defines the time width of the elementary detection in the time / frequency domain, also named pixels. A low value for TimeStep allows to monitor in greater detail the time variation of the wave parameters without interfering with the detection performance. Reducing TimeStep by a factor of 2 proportionately increases the number of pixels per family. In order to maintain a similar level of detection using TimeStep/2, ThresholdLFamMin must be multiplied by two and the standard deviations Sigma_t, Sigma_f, Sigma_v and Sigma_a can be reduced (see below).
  3. ThresholdConsistency (s)is the maximum consistency threshold for detection. This parameter is the first criterion to reject a pixel, and the choice of its value has a major impact on the detection performance. ThresholdConsistency has to be higher than the sampling period, otherwise the closure relation will never be fullfilled and no detection will be obtained. In unfavorable background noise conditions, including the consequent interferences in the inter-correlation functions, ThresholdConsistency must be raised. In addition, too high a value will noticeably increase the number of false detections linked to the incoherent nature of the background noise. The tuning of this value is achieved taking into account the frequency band being studied and the mean noise level prevailing at the station
  4. Q defines the threshold distance for integration of a sensor in the current subnetwork. There is integration if the distance between the sensor and the barycenter of the subnetwork is lower than Q. This threshold is useful to remove far-away sensor for which poorly-correlated signals are expected. If the maximum distance between sensors is not too large compared to the wavelength, Q can be set to a high value (typically 50). In that case, all sensors will be used for the calculation.
  5. ThresholdNbSensors is the minimum number of sensors participating to a detection. It ranges between 3 and the number of elements in the array. A high level of detection is obtained with a large value of ThresholdNbSensors.

2.3Frequency parameters