Multiresolution Analysis of a Seismological Signal

1. Filtering by Denoising

Fig. 1: original data.

The noise is clearly non-stationary - it is low up to the event near time point 1000, and then from time point 1500 onwards is is moderately low. We used a non-stationary noise model, to clip noise on each level of a multiscale transform. The transform used was the so-called à trous ("with holes") wavelet transform. The name à trous is due to the gaps used in the algorithm (which we won't go further into, here). This wavelet transform is unlike widely used transforms due to Mallat, Daubechies, Haar and others in that it is redundant: the user gets a set of versions of the data, each of which are the same length as the original input data. These are the sequence of resolution scales. The values at each resolution scale are termed wavelet coefficients. A final smooth version of the data (which we have termed residual, although one must be careful that this word is used in statistics - in regression analysis - in a very different sense) gives the following nice property: at each time point, we can add together the wavelet coefficient values at the successive resolution scales, plus the corresponding value in the smooth residual, to give exactly the input data value. That suggests immediately one way to filter the data: just drop one or more resolution scale value entirely, and recreate the data.

We could hope for a more sophisticated, statistically-based, alternative. This is to put small-valued wavelet coefficients to zero, at all resolution levels, if we feel that they are sufficiently small to be considered as noise only. That raises the question as to what noise is. Questions which this gives rise to include the following. How does noise arise in the first place? What is meant by a statistical noise model? How do we decide (statistically) what our guidelines are, in order to decide what is and what is not noise? What is the justification for carrying out denoising at the multiple resolution levels provided by a wavelet transform? An example of doing such denoising is shown in Fig. 2. This is not too bad in the regions to left and right of the earthquake event. But we haven't made much inroad into this event interval. The difference between original and filtered data (Fig. 3) shows this. In fact, our program has "protected" the earthquake event region, thinking that it comprised some very significant values. In a way this is true, and if we don't know how otherwise to filter the earthquake event interval, then it is maybe best to leave it alone.

Fig. 2: noise filtered data, based on a non-stationary noise model.

Fig. 3: the noise filtered data subtracted from the original data, indicating exactly what we have found as noise.

2. Wavelet Transform Visualization

Fig. 4: wavelet transform (using a redundant method, the B₃ spline à trous method).

Fig. 4 shows the succession of wavelet resolution scales (from upper left, to lower right) with the final smooth residual in the lower right. The number of scales (5 + residual here) is set by the user. An interesting thing happens by scale 5: the earthquake event in terms of time-interval is no longer significant enough to survive. It has almost disappeared.

To demarcate the earthquake event, we are looking for a sharp transition from what took place before, and maybe a sharp transition to normalcy which takes place afterwards. A term commonly applied to this in image and signal processing is edge detection. We are looking for an edge, a sharp transition. In statistical time series analysis, another term is sometimes used, a change-point. That term is more likely to imply a change from one state to another stable state. Our search for an edge, or multiple edges, in Fig. 1 is almost impossible, due to the noise in the data. Noise filtering, as seen in Fig. 2, did help somewhat, but not maybe as much as we would have liked.

We will find the different regime associated with the earthquake event by an unusual approach, which we have found to work well in similar circumstances. [It has also been used in (ref.) - I'm still looking...] If information exists on a number of resolution scale, as will be the case with very sharp exceptional information, then multiplying the information provided by these scales will accentuate the contrast. The multiplication is pixel-wise, or timestep-wise. The result of doing this for all scales (therefore no "selective ignorance" by removing scales, such as the smooth residual, which we could have done) is shown in Fig. 5. The ordinate (vertical axis of the plot) scale is huge. We are interested in the relativities, and these show up the region of interest. Also shown are indications of component shocks of the earthquake event, and after-shocks.

Fig. 5: the product of all scales, to accentuate edge-related contrasts.

To help the interpretation of this contrast-rich result in Fig. 5, Fig. 6 shows both it and the original data.

Fig. 6: the original data, Fig. 1, and the product result of Fig. 5, shown together.

Appendix

Steps followed: Given multiple column data, we used the programming language awk to quickly get just one column of values. The time steps corresponding to the values were not constant. We ignored this and took the data as if it had constant time increments. Then we translated it into FITS (an image or signal format, especially used in astronomy) using

 
im1d_convert 1b.dat 1b.fits 

 
s = readfits('1b.fits') 

Next we use the filtering program

 
mr1d_filter

 
mr1d_filter, s, sfilt, opt='-m 5'

s2000 = fltarr(2048)
s2000 = s + 2000.0
mr1d_filter, s2000, s2000filt, opt='-m 5'
s2000filt = s2000filt - 2000.0
plot,s2000filt
plot,s2000filt - s

mr1d_trans, s, strans, OPT="-n 6"
!p.multi=[0,3,2]
plot,strans.coef(*,0)
plot,strans.coef(*,1)
plot,strans.coef(*,2)
plot,strans.coef(*,3)
plot,strans.coef(*,4)
plot,strans.coef(*,5)

s1 = strans.coef(*,0)
s2 = strans.coef(*,1)
s3 = strans.coef(*,2)
s4 = strans.coef(*,3)
s5 = strans.coef(*,4)
s6 = strans.coef(*,5)
plot,s6*s5*s4*s3*s2*s1