Denoising ============================================= The noise in DAS signals mainly includes spike noise, common mode noise, stochastic noise and coherent noise. DASPy's data denoising module provides three functions: :ref:`Spike Removal`, :ref:`Common Mode Noise Removal` and :ref:`Stochastic Noise Removal`. Coherent noise can be removed by methods provided in :doc:`Wavefield Decomposition`. .. _Spike Removal: Spike Removal ------------------------------ Spikes are unusually large amplitudes and could be caused by laser frequency drift or laser noise. The spike removal function first applies the across-channel median filter and then the across-time median filter to generate a median map from the absolute amplitudes. Points with amplitudes exceeding a predefined threshold of the median map are identified as spikes. All spikes are subsequently substituted with interpolated values from adjacent channels. .. note:: The example data is the seismic signal recorded by Stanford DAS-1, which can be downloaded from ``_ , read and preprocess via: >>> from daspy import read, DASDateTime >>> sec = read('Stanford1_20180104_113058.566+0000.sgy', ch2=320) UserWarning: This data format segy doesn't include channel interval. Please set Section.dx manually. >>> sec.dx = 8 >>> sec.start_time = DASDateTime(2018, 1, 4, 11, 30, 58, 566000) >>> origin_time = DASDateTime(2018, 1, 4, 11, 34, 44) >>> sec.bandpass(1, 20) >>> sec.trimming(tmin=origin_time+2, tmax=origin_time+12) :: >>> sec.plot(xmode='channel') .. image:: ../media/spike_removal0.png :width: 600 :: >>> sec.spike_removal() >>> sec.plot(xmode='channel') .. image:: ../media/spike_removal1.png :width: 600 .. _Common Mode Noise Removal: Common Mode Noise Removal ------------------------------ Common-mode noise, also known as in-phase noise is generated by vibrations of the optoelectronic system and arises on all channels simultaneously. DASPy employs spatial median or mean of waveforms to obtain common mode noise. Subsequently, we compute the correlation coefficient with the channel record and the common-mode noise, multiply the common-mode noise by the coefficient, and subtract it from the channel record. .. note:: The example data is the waveforms recorded by channels away from the coast of the RAPID dataset, which can be downloaded from ``_ , read and preprocess via: >>> import numpy as np >>> from daspy import read >>> sec = read('North-C1-LR-P1kHz-GL50m-Sp2m-FS200Hz_2021-11-04T015902Z.h5') >>> sec.trimming(mode=0, xmin=18000, xmax=30000) >>> sec.scale = 2 * np.pi / 2 ** 16 # data scaling factor, see sec.headers['Acquisition']['Raw[0]']['attrs']['RawDataUnit'] >>> sec.phase2strain(1550.12 * 1e-9, 0.78, sec.headers['Acquisition']['Custom']['attrs']['Fibre Refractive Index']) # convert optical phase shift to strain >>> sec.bandpass(15, 27, detrend=True, taper=0.1) >>> sec.trimming(tmin=sec.start_time+20, tmax=sec.start_time+30) :: >>> sec.plot() .. image:: ../media/common_mode_noise_removal0.png :width: 600 :: >>> sec.common_mode_noise_removal() >>> sec.plot() .. image:: ../media/common_mode_noise_removal1.png :width: 600 .. _Stochastic Noise Removal: Stochastic Noise Removal ------------------------------ The inherent stochastic noise in DAS data is primarily caused by instrumental deficiencies such as sampling error and phase noise. DASPy remove stochastic noise by curvelet transform. .. note:: Consistent with the example data used by :ref:`Spike Removal` . The waveform with spikes removed is used here. Read and preprocess via: >>> from daspy import read, DASDateTime >>> sec = read('Stanford1_20180104_113058.566+0000.sgy', ch2=320) UserWarning: This data format segy doesn't include channel interval. Please set Section.dx manually. >>> sec.dx = 8 >>> sec.start_time = DASDateTime(2018, 1, 4, 11, 30, 58, 566000) >>> sec.bandpass(1, 20) >>> origin_time = DASDateTime(2018, 1, 4, 11, 34, 44) >>> sec.trimming(tmin=origin_time-10, tmax=origin_time+12) >>> sec.spike_removal() Use a noise record as the noise baseline and remove the noise with a soft threshold (default) in the curvelet domain: >>> sec_eq = sec.copy().trimming(tmin=origin_time+2, tmax=origin_time+12) # earthquake records >>> sec_ns = sec.copy().trimming(tmin=origin_time-10, tmax=origin_time) # noise records >>> sec_eq_soft = sec_eq.copy().curvelet_denoising(noise=sec_ns) Also using the reference noise record, removing the noise with a hard threshold in the curvelet domain can keep the absolute amplitude of the waveform unchanged and cause little distortion: >>> sec_eq_hard = sec_eq.copy().curvelet_denoising(noise=sec_ns, soft_thresh=False) When there is no reference noise record available, the function will calculate the inflection points of the curvelet coefficients to determine the noise threshold. It is recommended to set ``pad=0`` and adjust the ``knee_fac`` parameter to reduce artificial artifacts (this method is not recommended) : >>> sec_eq_knee = sec_eq.copy().curvelet_denoising(pad=0, knee_fac=0.1) Plot the original waveform and the above three denoising effects: >>> import matplotlib.pyplot as plt >>> fig, ax = plt.subplots(2, 2, figsize=(6,6), sharex=True, sharey=True, dpi=200) >>> sec_eq.plot(ax=ax[0,0], xmode='channel', vmax=0.2, xlabel=False, colorbar=False) >>> sec_eq_soft.plot(ax=ax[0,1], xmode='channel', vmax=0.2, xlabel=False, ylabel=False, colorbar=False, title='soft thresh') >>> sec_eq_hard.plot(ax=ax[1,0], xmode='channel', vmax=0.2, colorbar=False, title='hard thresh') >>> sec_eq_knee.plot(ax=ax[1,1], xmode='channel', vmax=0.2, ylabel=False, colorbar=False, title='without reference noise') >>> plt.tight_layout() >>> plt.show() .. image:: ../media/curvelet_denoising.png :width: 700