Home | About | Journals | Submit | Contact Us | Français |

**|**Scientific Reports**|**PMC5607239

Formats

Article sections

- Abstract
- Introduction
- Morphological Decomposition
- Traditional morphological reconstruction
- Orthogonalized morphological reconstruction
- Solution of orthogonalization operator
- Implementation of the orthogonalized morphological reconstruction
- Efficiency and effectiveness analysis of the orthogonalized morphological reconstruction
- Test of the orthogonalized morphological reconstruction
- Application to a real data set
- Conclusion
- Electronic supplementary material
- References

Authors

Related links

Sci Rep. 2017; 7: 11996.

Published online 2017 September 20. doi: 10.1038/s41598-017-09711-2

PMCID: PMC5607239

Yangkang Chen, Email: moc.liamg@6102kynehc.

Received 2017 April 6; Accepted 2017 July 28.

Copyright © The Author(s) 2017

Microseismic method is an essential technique for monitoring the dynamic status of hydraulic fracturing during the development of unconventional reservoirs. However, one of the challenges in microseismic monitoring is that those seismic signals generated from micro seismicity have extremely low amplitude. We develop a methodology to unveil the signals that are smeared in the strong ambient noise and thus facilitate a more accurate arrival-time picking that will ultimately improve the localization accuracy. In the proposed technique, we decompose the recorded data into several morphological multi-scale components. In order to unveil weak signal, we propose an orthogonalization operator which acts as a time-varying weighting in the morphological reconstruction. The orthogonalization operator is obtained using an inversion process. This orthogonalized morphological reconstruction can be interpreted as a projection of the higher-dimensional vector. We first test the proposed technique using a synthetic dataset. Then the proposed technique is applied to a field dataset recorded in a project in China, in which the signals induced from hydraulic fracturing are recorded by twelve three-component (3-C) geophones in a monitoring well. The result demonstrates that the orthogonalized morphological reconstruction can make the extremely weak microseismic signals detectable.

It has been shown that microseismic monitoring has a significant potential to characterize physical processes related to fluid injections and extractions in hydrocarbon and geothermal reservoirs^{1,2}. In general the microseismicity is recorded by downhole or shallow surface geophone arrays, which offers the significant advantages of being sufficiently close to the fracture and being unaffected by the free surface^{3}. There are two main physical processes involved in hydraulic fracturing: 1) penetration of the injected fluid into the pre-existing cracks and pore spaces when the injection pressure is lower than the minimum compressive stress, and 2) opening of new fractures when the injection pressure is high enough. The events generated during injection and also after injection can occur over hours^{4}. Localization of the associated microseismic events enables imaging of the fracture network. This technique has been widely studied and applied in petroleum and gas exploration^{1,5–11}, and mining engineering^{12–15}. However, an inevitable problem existing in the microseismic monitoring is that the energy stimulated from the hydraulic fracturing is extremely weak, compared with the background noise^{16}. The weak signal is easily masked, resulting in loss of microseismic events. A poor signal-to-noise ratio (S/N) can lead to unauthentic arrival time-picks^{17} and localization of microseismic events^{18}. All of these will negatively affect the performance of microseismic monitoring and resulted fracture imaging^{19}, as well as solving source mechanisms^{20}. Improving the S/N will ultimately improve the microseismic event detection. In microseismic monitoring, the most commonly used method for attenuating background noise and detecting weak signal is frequency filtering^{21}. However, frequency filtering typically fails in separating noise and signal when they share the same frequency band. Researchers put a lot of effort into the noise suppression problem^{22}, and developed different techniques using different approaches such as: median filtering^{23}, various kinds of mathematical transform based approaches^{24–26}, and matrix completion based approaches^{27,28}. In addition, Kong *et al*. develop ed a nonlinear signal detector, which passes only signals showing spatial coherence and having slowness within an allowed range. Schimmel and Paulssen^{30} use d an instantaneous phase based amplitude-unbiased coherency measure, weighting the samples of an ordinary, linear stack, to detect weak signals in global seismology. Gibbons and Ringdal^{31} illustrated the power of an array-based waveform correlation approach by detecting the low magnitude seismic events in the 1997 August 16 Kara Sea event. Mousavi *et al*.^{32} proposed a simultaneous microseismic denoising and onset detection technique based on the synchrosqueezed continuous wavelet transform and custom thresholding of single-channel data. Mousavi and Langston^{33} designed fast algorithm for noise level estimation and noise reduction of micro-seismic data, using minimally controlled recursive averaging and neighborhood shrinkage estimators.

The denoising approach of this paper is based on mathematical morphological decomposition and reconstruction^{34,35}. Mathematical morphology is a nonlinear methodology for the analysis and processing of geometrical structures. Matheron^{34} describe d the random set integral geometry theory and topological logic theories thoroughly and set up a consistent foundation for mathematical morphology. Later on, Serra^{35} suggested the theory and method of mathematical morphology which was widely applied in two-value image processing. Then, Koskinen *et al*.^{36} introduce d the soft mathematical operations, which can maintain most of the properties of standard morphological operations. Sinha and Dougherty^{37} developed a generalization of binary mathematical morphology based on fuzzy set theory, in which images are modeled as fuzzy subsets of the Euclidean plane or Cartesian grid, and the morphological operations are defined in terms of a fuzzy index function. The mathematical morphological filtering (MMF) was first introduced into seismic data processing by Wang *et al*.^{38}. Unlike traditional methods in seismic data processing, the basis of MMF are logical operation and set theory, which can provide us a tool to process signal over the complete frequency bandwidth, improving the S/N and maintaining the resolution. Later, this method rapidly developed and was widely applied in seismic data processing. For example, Li *et al*.^{21} proposed a compound top-hat filter (CTF) extracting the large-scale information by combining opening and closing operations, and subsequently subtracting it from the microseismic data.

In this paper, we further develop a seismic application of the mathematical morphology and propose a multi-scale morphological decomposition based method to unveil weak signal in microseismic monitoring. In order to unveil weak signal, an orthogonalization operator is proposed and introduced into the process of multi-scale morphological reconstruction. The mathematical nature of the proposed orthogonalization operator is a projection operator that projects the input signal on a sub-space spanned by several selected morphological basis vectors. The assumption for this approach is that the weak signal is orthogonal to the background noise. Unlike the traditional morphological reconstruction approaches, the orthogonalized morphological reconstruction transforms the reconstruction problem into an inversion problem. However, like most of the inversion problems in geophysics, this inversion problem is ill-posed. A regularization (or a penalty) term is necessary to optimally stabilize the objective function. In this study we use the shaping regularization approach^{39} that is more convenient for solving the inversion problem in orthogonalized morphological reconstruction compared to other regularization technique such as Tikhonov’s method^{40}. In the following sections, after explaining the methodology we first conducts synthetic data experiments to test the performance of the proposed orthogonalized morphological reconstruction approach. Then the proposed approach is applied to a real 3-C microseismic data set. Compared with the state-of-the-art algorithms, the proposed approach demonstrates a superior performance.

Mathematical morphology^{34,35}, is a well known nonlinear image processing method, which was originally motivated from the research of the relation between the penetrability of a porous medium and its lamination. It starts as a set theoretical approach for the analysis of geometrical structures but can also deal with both function and set in the Euclidean space. A morphological operation is the interaction of an objective set or function with another set or function called structuring element (SE). The morphological scale of the SE determines the scale information of the signal that is extracted under such an operation. The morphological scale can be conceptually understood as the relative structure. Let **d** be a seismic trace and **f** ⊆ *ℝ* be a set of amplitude values. The value of a sample *t* in **d** is represented by *d*(*t*) ∈ **f**. The morphological dilation *φ*_{b}(**d**) and erosion *φ*_{b}(**d**) are the morphological operations that process **d** with the SE *b*(*τ*) as^{41}:

$${\phi}_{b}(\mathbf{d})=\underset{\tau}{\vee}b(\tau )+d(t-\tau ),$$

1

$${\varphi}_{b}(\mathbf{d})=\underset{\tau}{\wedge}b(\tau )-d(\tau -t),$$

2

where ∨ denotes supremum, and ∧ denotes infimum. Both *t* and *τ* are samples. It can be seen that the morphological dilation is an operation that “grows” or “thickens” the object, while the morphological erosion is an operation that “shrinks” or “thins” the object. The sequential combination of the morphological erosion (or dilation) and morphological dilation (or erosion) creates the morphological opening *χ*_{b}(**d**) (or closing *ψ*_{b}(**d**)) as:

3

$${\psi}_{b}(\mathbf{d})={\varphi}_{b}({\phi}_{b}(\mathbf{d}\mathrm{)).}$$

4

We now use morphological opening and closing to represent data **d**. Consider {*χ*_{bk}(**d**)}, $k\in \mathrm{[1,}\phantom{\rule{.25em}{0ex}}K]$ and {*ψ*_{bk}(**d**)}, $k\in \mathrm{[1,}\phantom{\rule{.25em}{0ex}}K]$, two indexed families of morphological opening and closing, respectively. Typically, the index *k* denotes the morphological scale. Whereupon, **d** can be represented as:

$$\mathbf{d}=\sum _{k\mathrm{=1}}^{K}[{\chi}_{{b}_{k-1}}({\psi}_{{b}_{k-1}}(\mathbf{d}))-{\chi}_{{b}_{k}}({\psi}_{{b}_{k}}(\mathbf{d}))]+{\chi}_{{b}_{K}}({\psi}_{{b}_{K}}(\mathbf{d})),$$

5

or

$$\mathbf{d}=\sum _{k\mathrm{=1}}^{K}[{\psi}_{{b}_{k-1}}({\chi}_{{b}_{k-1}}(\mathbf{d}))-{\psi}_{{b}_{k}}({\chi}_{{b}_{k}}(\mathbf{d}))]+{\psi}_{{b}_{K}}({\chi}_{{b}_{K}}(\mathbf{d}\mathrm{)).}$$

6

Equations (^{5}) and (^{6}) can also be written as:

$$\begin{array}{ll}\mathbf{d}\hfill & =\frac{1}{2}(\mathbf{d}+\mathbf{d})\hfill \\ \hfill & =\frac{1}{2}\{{\psi}_{{b}_{K}}({\chi}_{{b}_{K}}(\mathbf{d}))+{\chi}_{{b}_{K}}({\psi}_{{b}_{K}}(\mathbf{d}))\hfill \\ \hfill & \phantom{\rule{1em}{0ex}}+\sum _{k\mathrm{=1}}^{K}[{\chi}_{{b}_{k-1}}({\psi}_{{b}_{k-1}}(\mathbf{d}))-{\chi}_{{b}_{k}}({\psi}_{{b}_{k}}(\mathbf{d}))+{\psi}_{{b}_{k-1}}({\chi}_{{b}_{k-1}}(\mathbf{d}))-{\psi}_{{b}_{k}}({\chi}_{{b}_{k}}(\mathbf{d}\mathrm{))]\}.}\hfill \end{array}$$

7

So far, the initial data **d** is represented by an additive decomposition with *K*+1 scales. Figure 1 gives an example of the morphological decomposition of a Ricker wavelet with 7 scales. The 1st trace (scale 0) is the initial wavelet. The 2*nd*–8th traces are the 7 scale components.

For convenience, let:

$${\mathbf{c}}_{k}=\{\begin{array}{ll}\frac{1}{2}\{{\chi}_{{b}_{k-1}}({\psi}_{{b}_{k-1}}(\mathbf{d}))-{\chi}_{{b}_{k}}({\psi}_{{b}_{k}}(\mathbf{d}))+{\psi}_{{b}_{k-1}}({\chi}_{{b}_{k-1}}(\mathbf{d}))-{\psi}_{{b}_{k}}({\chi}_{{b}_{k}}(\mathbf{d}))]\},\hfill & \phantom{\rule{.1em}{0ex}}k\in \phantom{\rule{.5em}{0ex}}\mathrm{[1,}\phantom{\rule{.25em}{0ex}}K],\phantom{\rule{.1em}{0ex}}\hfill \\ \frac{1}{2}\{{\psi}_{{b}_{K}}({\chi}_{{b}_{K}}(\mathbf{d}))+{\chi}_{{b}_{K}}({\psi}_{{b}_{K}}(\mathbf{d}))\},\hfill & \phantom{\rule{.1em}{0ex}}k=K+\mathrm{1,}\phantom{\rule{.1em}{0ex}}\hfill \end{array}$$

8

where **c**_{k}, $k\in \mathrm{[1,}\phantom{\rule{.25em}{0ex}}K+\mathrm{1]}$, are the morphological multi-scale components. The value of a sample *t* in **c**_{k} is represented by *c*_{k}(*t*) ∈ **f**. The nature of the multi-scale morphological decomposition is to decompose the discrete data set **d** into a series of primary subsets **c**_{k}, which satisfies that:

$$\mathbf{d}=\underset{k\mathrm{=1}}{\overset{K+1}{\cup}}{\mathbf{c}}_{k},$$

9

10

where denotes empty set. The reconstruction of data by **c**_{k} can be represented as:

$$\sum _{{\mathbf{c}}_{k}\in \mathbf{E}}[{\mathbf{\sigma}}_{k}{\mathbf{c}}_{k}]=\mathbf{d},$$

11

where **E** is a subset of ${\{{\mathbf{c}}_{k}\}}_{k\in \mathrm{[1,}K]}$, that $\mathbf{E}\underset{\_}{\underset{\_}{\subset}}{\{{\mathbf{c}}_{k}\}}_{k\in \mathrm{[1,}K]}$. Constant ${\mathrm{\sigma}}_{k}\in \mathrm{[0,}\phantom{\rule{.25em}{0ex}}\mathrm{1]}$ is the weighting coefficient that controls energy from different scale components. This decomposition allows for full reconstruction of the original data, when $\mathbf{E}=\{{\mathbf{c}}_{k}{\}}_{k\in \mathrm{[1,}K]}$ and σ_{k} ≡ 1.

In traditional morphological reconstruction, the weighting coefficient σ_{k} is chosen manually, which makes the reconstruction subjective. In addition, it is difficult to choose an appropriate weighting coefficient, and the choosing process costs a lot of time and manual endeavor. For weak signal detection, we define the orthogonalized morphological reconstruction, by changing σ_{k} from simple constant to a more flexible operator, in other words, allowing σ_{k} to change with *t*:

$$\mathrm{arg}\underset{{\sigma}_{k}(t)}{\mathrm{min}}\sum _{{\mathbf{c}}_{k}\in \mathbf{E}}{\Vert {\mathbf{\sigma}}_{k}(t){c}_{k}(t)-d(t)\Vert}_{F}^{2},$$

12

where $\parallel \phantom{\rule{-.25em}{0ex}}\cdot \phantom{\rule{-.25em}{0ex}}{\parallel}_{F}^{2}$ represents the squared Frobenius norm of a function. Equation (^{12}) can be also represented as:

$$\mathrm{arg}\underset{{\mathbf{\Sigma}}_{k}}{\mathrm{min}}\sum _{{\mathbf{c}}_{k}\in \mathbf{E}}{\Vert {\mathbf{\Sigma}}_{k}{\mathbf{c}}_{k}-\mathbf{d}\Vert}_{F}^{2},$$

13

where **Σ**_{k} is a diagonal matrix composed by **σ**_{k}(*t*): **Σ**_{k} = *d**i**a**g*(**σ**_{k}(*t*)). Thus, the orthogonalized morphological reconstruction holds as:

$$\sum _{{\mathrm{c}}_{k}\in \mathbf{E}}[{\mathbf{\Sigma}}_{k}{\mathbf{c}}_{k}]=d\mathrm{.}$$

14

The geometrical nature behind equation (^{13}) is a projection of the higher-dimensional vector, i.e., the initial data **d**, on a lower-dimensional space spanned by several selected morphological basis vectors. Figure 2 gives a diagrammatic drawing. Vector $\overrightarrow{\mathbf{c}}$ is on the line *l*. Operation **Σ** is a stretching transformation acting on vector $\overrightarrow{\mathbf{c}}$. Equation (^{13}) is actually to find the projection of vector $\overrightarrow{\mathbf{d}}$ on line *l* in the least-squares sense. Thus, an orthogonal decomposition of $\overrightarrow{\mathbf{d}}$ holds as:

$$\overrightarrow{\mathbf{d}}=\mathbf{\Sigma}\overrightarrow{\mathbf{c}}+(\overrightarrow{\mathbf{d}}-\mathbf{\Sigma}\overrightarrow{\mathbf{c}}),$$

15

$$\overrightarrow{\mathbf{0}}=\mathbf{\Sigma}\overrightarrow{\mathbf{c}}\cdot (\overrightarrow{\mathbf{d}}-\mathbf{\Sigma}\overrightarrow{\mathbf{c}}),$$

16

where ⋅ denotes Hadamard (or Schur) product. Hence, we name **Σ**
*orthogonalization operator*. If we consider $\mathbf{\Sigma}\overrightarrow{\mathbf{c}}$ as signal $\overrightarrow{\mathbf{s}}$, and accordingly $\overrightarrow{\mathbf{d}}-\mathbf{\Sigma}\overrightarrow{\mathbf{c}}$ as background noise $\overrightarrow{\mathbf{n}}$, equations (^{15}) and (^{16}) become the classical models used in^{42–45},

$$\overrightarrow{\mathrm{d}}=\overrightarrow{\mathbf{s}}+\overrightarrow{\mathbf{n}},$$

17

$$\overrightarrow{\mathbf{0}}=\overrightarrow{\mathbf{s}}\cdot \overrightarrow{\mathbf{n}}\mathrm{.}$$

18

Therefore, if we assume the weak signal is orthogonal to the background noise in microseismic monitoring, $\mathbf{\Sigma}\overrightarrow{\mathbf{c}}$ is an estimation of the weak signal.

The inversion problem in equation (^{13}), however, is ill-posed. To stabilize the optimization, an extra regularization term is necessary to solve equation (^{13}):

$$\mathrm{arg}\underset{{\mathbf{\Sigma}}_{k}}{\mathrm{min}}\sum _{{\mathbf{c}}_{k}\in \mathbf{E}}[{\Vert {\mathbf{\Sigma}}_{k}{\mathbf{c}}_{k}-\mathbf{d}\Vert}_{F}^{2}+\mathcal{R}({\mathbf{\Sigma}}_{k})],$$

19

where *ℛ* represents the regularization operator. For convenience, we rewrite equation (^{19}) as:

$$\mathrm{arg}\underset{{\mathbf{\sigma}}_{k}}{\mathrm{min}}\sum _{{\mathbf{c}}_{k}\in \mathbf{E}}[{\Vert {\mathbf{C}}_{k}{\mathbf{\sigma}}_{k}-\mathbf{d}\Vert}_{F}^{2}+\mathcal{R}({\mathbf{\sigma}}_{k})],$$

20

where **C**_{k} is a diagonal matrix composed by *c*_{k}(*t*): **C**_{k} = *d**i**a**g*(*c*_{k}(*t*)). **σ**_{k} is a column vector composed by **σ**_{k}(*t*): **σ**_{k} = [**σ**_{k}(*t*)]^{T}. Note that **Σ**_{k}*c*_{k} = **C**_{k}**σ**_{k} = **σ**_{k} ⋅ **c**_{k}.

One of the most commonly used regularization approaches is Tikhonov’s regularization^{40}, in which one additionally attempts to minimize the norm of **T***σ*_{k}, where **T** is the regularization operator^{39}. The regularized problem can be expressed as:

$$\mathrm{arg}\underset{{\mathbf{\sigma}}_{k}}{\mathrm{min}}\sum _{{\mathbf{C}}_{k}\in \mathbf{E}}[{\Vert {\mathbf{C}}_{k}{\mathbf{\sigma}}_{k}-\mathbf{d}\Vert}_{F}^{2}+{\epsilon}^{2}{\Vert \mathbf{T}{\mathbf{\sigma}}_{k}\Vert}_{F}^{2}],$$

21

where *ε* is a scalar scaling parameter. The formal solution has the well-known form,

$${\overline{\mathbf{\sigma}}}_{k}={({\mathbf{C}}_{k}^{T}{\mathbf{C}}_{k}+{\epsilon}^{2}{\mathbf{T}}^{T}\mathbf{T})}^{-1}{\mathbf{C}}_{k}^{T}\mathbf{d}\mathrm{.}$$

22

where ${\overline{\mathbf{\sigma}}}_{k}$ present a least-squares estimate of *σ*_{k}. Tikhonov’s regularization can be interpreted as a roughening approach. Although Tikhonov’s regularization is effective, the parameter *ε* and regularization operator **T** are typically difficult to choose, from the user perspectives^{46,39} proposed a particularly convenient shaping regularization approach, in which, a triangle shaping operator **Γ** is proposed and introduced into the iterative inversion as a fundamental operation. The relation between regularization operator **T** and shaping operator **Γ** can be expressed as^{39}:

23

Combining equation (^{22}) and (^{23}), we have:

$${\overline{\mathbf{\sigma}}}_{k}={[\mathbf{I}+\mathbf{\Gamma}({\mathbf{C}}_{k}^{T}{\mathbf{C}}_{k}-\mathbf{I})]}^{-1}\mathbf{\Gamma}{\mathbf{C}}_{k}^{T}\mathbf{d}\mathrm{.}$$

24

By introducing scaling of A by $\mathrm{1/}\lambda $ in equation (^{24}), we can rewrite it as:

$${\overline{\mathbf{\sigma}}}_{k}={[{\lambda}^{2}\mathbf{I}+\mathbf{\Gamma}({\mathbf{C}}_{k}^{T}{\mathbf{C}}_{k}-{\lambda}^{2}\mathbf{I})]}^{-1}\mathbf{\Gamma}{\mathbf{C}}_{k}^{T}\mathbf{d}\mathrm{.}$$

25

where *λ* is an introduced parameter controlling the physical dimensionality and enabling fast convergence when inversion is implemented iteratively^{27}.

The SE plays an important role in the morphological decomposition and reconstruction. The SE has three parameters: shape, height (the amplitude of SE), and width (the width of definitional domain of SE). Generally speaking, the shape of SE can be a semicircle, a triangle, or a straightline. The SEs with different parameters has different scales. When the shape of a SE is fixed, its scale increases as the height decreases (or as the width increases). A SE with a large (or small) scale indicates that it has a fat (or slim) structure (i.e., its shape is close to the shape of a constant (or *δ*) function). The comparison of scale among the three shapes is as follow:

$$Scale(straightline)>Scale(semicircle)>Scale(triangle\mathrm{).}$$

26

In the morphological decomposition, we need a series of SE with different scales to obtain the different morphological information of the input data. For a specific morphological decomposition, a commonly used strategy to produce the SE family **b**_{K} is that we fix the shape of the SE and gradually increase both its height and width to produce different SEs. The rate of increase determines the performance of decomposition. Another more convenient strategy is that the *i* th SE *b*_{i} can be produced by *i* − 1 times self morphological dilation:

$${\mathbf{b}}_{i}=\underset{i-1\phantom{\rule{.5em}{0ex}}times}{\underset{\u23df}{{\phi}_{{b}_{1}}\mathrm{(...}{\phi}_{{b}_{1}}({\phi}_{{b}_{1}}}}({\mathbf{b}}_{1}\mathrm{))).}$$

27

An iterative optimization can greatly improve efficiency in solving an inverse problem when the computational scale is large. We choose the classical conjugate gradient method^{47} to iteratively implement the orthogonalized morphological reconstruction approach. The conjugate gradient algorithm requires symmetric positive definite operators. So the shaping operator splits into two matrices, **Γ** = **H****H**^{T}. Equation (^{25}) can then be written as:

$${\overline{\mathbf{\sigma}}}_{k}=\mathbf{H}{[{\lambda}^{2}\mathbf{I}+{\mathbf{H}}^{T}({\mathbf{C}}_{k}^{T}{\mathbf{C}}_{k}-{\lambda}^{2}\mathbf{I})\mathbf{H}]}^{-1}{\mathbf{H}}^{T}{\mathbf{C}}_{k}^{T}\mathbf{d}\mathrm{.}$$

28

The estimated weak signal **s** by the orthogonalized morphological reconstruction can be represented as:

$$\mathbf{s}\approx \sum _{{\mathbf{c}}_{k}\in \mathbf{E}}[{\mathbf{C}}_{k}{\overline{\mathbf{\sigma}}}_{k}\mathrm{].}$$

29

The proposed technique first decomposes the input data into a series of components with different morphological features, and then reconstructs the signal by several selected components with an orthogonalization operator. Decomposition with a higher order can obtain a more careful multi-scale morphology analysis of the input data, and accordingly is easier to separate signal and noise. Unfortunately, a large number of decompositions will pose a very expensive computational cost. Our experience shows that 4–10 decomposed components are appropriate for most seismic data sets, taking the compromise between efficiency and effectiveness into consideration. Figure 3 demonstrates an experimental analysis of the proposed orthogonalized morphological reconstruction method. Figure 3(a,b) show the computing time costs and denoising performance analysis varying with different numbers of decomposition of the input data. We can observe that, as the decomposition number increases, the computational time increases. The denoising performance of the proposed technique is reinforced as the decomposition number increases within a relatively small value (2–6), but maintains relatively stable when the decomposition number is greater than 6. Figure 3(c) shows the denoising performance varying with different input S/Ns.

A synthetic signal is used to test our proposed method in this section. The first experiment is shown in Fig. 4. The synthetic signal is a Ricker wavelet with 100 Hz dominant frequency and $\pi \mathrm{/2}$ initial phase, as shown by the 1st trace. The synthetic noise is broadband Gaussian noise as shown in the 2nd trace. The 1st trace is added with the 2nd trace as the input data (the 3th trace). The S/N of the input data is −11.6971 dB and the definition of S/N is shown below^{48}:

$$S/N\phantom{\rule{.15em}{0ex}}=10\phantom{\rule{.15em}{0ex}}\mathrm{log}\phantom{\rule{.10em}{0ex}}\frac{{\Vert \mathbf{s}\Vert}_{F}^{2}}{{\Vert \mathbf{n}\Vert}_{F}^{2}},$$

30

where **s** is the true signal and **n** denotes the added noise. We can see that the signal is masked by the strong background noise and hardly to detect. For detection of the masked signal, the input data is decomposed into seven multi-scale morphological components as shown in the 4th–10th traces. It can be observed that most energy of noise is decomposed into the 1st and 2nd multi-scale components (the 4th and 5th traces) and the signal can be followed more or less in the rest components. Thus the 3th–7th multi-scale components (the 6th–10th traces) are used to reconstruct the signal. The results using the proposed and conventional morphological reconstruction approaches are shown in the 11th and 12th traces, respectively. The weighting coefficients in conventional approach are chosen manually as $\mathrm{(1,1,1,1,0.5)}$ associated to the 6th–10th traces taking the compromise between preserving signal and suppressing noise into consideration. It is clear that the signal is much more detectable after both two reconstructions. But the proposed approach suppress more background noise and makes the signal easier to detect than the conventional. As a comparison, the commonly used band-pass filtering is applied to the input data. The results using three different band-pass filterings (trapezoidal bands 10–20–180–190 Hz, 70–80–120–130 Hz and 60–70–170–180 Hz) are plotted in the 13th, 14th and 15th traces. The weak band-pass filtering can preserve signal well but pass a lot of background noise at the same time. The strong band-pass filtering can suppress more noise but damage the signal. The corresponding errors (differences between the true signal and processed results) associated to traces 11–15 are presented in traces 16–20 respectively. It is obvious that the proposed orthogonalized morphological reconstruction obtains the the smallest error. The S/Ns of the processed results using the proposed and conventional morphological reconstruction approaches and three band-pass filterings are 10.8905, −0.1852, −0.1788, 0.9289 and 0.6348 dB, respectively. The calculation of S/N refers to equation (^{30}), except that *n* denotes the error. It is obvious that the proposed approach obtain the highest S/N. The cross-correlation coefficients between original signal and the five denoised signals are 0.9585, 0.6446, 0.7086, 0.6163, and 0.6913, respectively.

The first synthetic example. From left to right: the 1st trace: signal, the 2nd trace: Gaussian noise, the 3th trace: signal + Gaussian noise (input data), the 4th–10th traces: seven multi-scale components, the 11th trace: orthogonalized morphological **...**

Figure 5 shows time-frequency spectra of clean data, noisy data, two reconstructions and three filtered results, which give us a more detailed comparison. The time-frequency spectrum is obtained by using standard Stockwell transform^{49}. As we can see from Fig. 5(a), the energy of the synthetic signal is concentrated in the area of 0.12–1.07 ms and 0–300 Hz. Contaminated by Gaussian noise (Fig. 5(b)), the time-frequency spectrum become noisy. The result using the proposed technique is shown in Fig. 5(c). An excellent reconstruction of the initial synthetic signal is obtained. Most of the noise has been attenuated and the signal’s energy is much more distinct. Figure 5(d) presents the reconstructed data obtained after using the conventional morphological reconstruction technique. Even though the quality of data is improved, the detected signal is still strongly affected by the noise. Figure 5(e,f), show the filtered data by different band-pass filters. The band-pass filtering is achieved by the low-cutoff and high-cutoff in frequency domain. As we can observe from Fig. 5(e,f), the low and high frequency components are removed but the mixed parts in the same frequency band still exist. The starting time of the detected signals in Fig. 5(d,e,f), are not as clear as that shown in Fig. 5(c), indicating that the time-picking would be better performed on the record processed by the proposed orthogonalized morphological reconstruction technique.

Comparison of time-frequency spectrums of (**a**) clean data, (**b**) noisy data, (**c**) orthogonalized morphological reconstruction, (**d**) conventional morphological reconstruction, (**e**) filtered data (10–20–180–190 Hz), (**f**) filtered data (70–80–120–130 **...**

The second example is demonstrated in Fig. 6. The synthetic signal (the 1st trace) is same to that in the first experiment. The added background noise consists of Gaussian noise (the 2nd trace) and limited band (40–160 Hz) random noise (the 3th trace). The input data is the sum of the 1st, 2nd and 3th traces as shown in the 4th trace. The S/N of the input data is −12.5386 dB. Similarly, the input data is decomposed into seven multi-scale morphological components as shown in the 5th–11th traces. In this experiment, we choose the 3 th–6th multi-scale components (the 7th–10th traces) to reconstruct the signal, taking the compromise between signal preservation and noise removal into consideration. The proposed and conventional reconstructions are plotted in the 12th and 13th traces. The weighting coefficients in conventional approach are chosen manually as $\mathrm{(1,1,1,1)}$ associated to the 7th–10th traces. It can be seen that, both approaches improve the detectability of the signal, but the proposed approach gives a better result. The three filtered data using band-pass filtering, with trapezoidal bands 10–20–180–190Hz, 70–80–120–130Hz and 60–70–170–180Hz, are shown in the 14th, 15th and 16th traces. The results are unacceptable. We still hardly detect the signal in the filtered data. The S/Ns of the processed results using proposed and conventional morphological reconstruction approaches and three band-pass filterings are 4.9067, −2.8560, −6.2078, −3.0327, and −5.6471 dB, respectively. The cross-correlation coefficients between original signal and the five denoised signals are 0.8254, 0.4392, 0.3944, 0.4095, and 0.3838, respectively. Similarly, the time-frequency spectr a of clean data, noisy data, two reconstructions and two filtered results are shown in Fig. 7. The manually added limited-band random noise increases the difficulty of detecting the weak signal. As we can observe from Fig. 7(b), the time-frequency spectrum is extremely noisy and particularly several energy clusters in the area of 0.17–0.3ms and 0–200Hz can seriously obstruct the detection of the true signal. Denoising by using the orthogonalized morphological reconstruction technique leads to the results depicted in Fig. 7(c). The noise is clearly suppressed. By comparing the true signal and reconstructed signal, we can see that the two signals are very similar except for a slight amplitude damage. However, the signal would be easier to pick than before, and the slight amplitude damage is not significant, considering the totally removed noise and the observable useful signal s. The conventional morphological reconstruction approach also remove some noise, but the noise energy clusters are still noticeable. Figure 7(e,f,g) demonstrate the three filtered results. As expected, by using band-pass filtering, noise that shares the same frequency band with signal cannot be separated.

The second synthetic example. From left to right: the 1st trace: signal, the 2nd trace: Gaussian noise, the 3th trace: limited band random noise, the 4th trace: signal + Gaussian noise + limited band random noise (input data), the 5th–11th **...**

The proposed orthogonalized morphological reconstruction is applied to a real microseismic monitoring dataset recorded in the west of China. There are twelve downhole 3-C geophones to monitor seismic activity. There are eight injection stages in this project. The magnitude of microseismic events ranges from −3.86 to −0.135 Mw. The data used in this study is produced in the last stage, in which the recording time is the longest in the whole project. Section “Supplementary material” gives the detailed information for this dataset. In this dataset, the signal induced from hydraulic fracturing is very weak when the signal reach the receivers. A lot of useful signal s cannot be detected immediately, which leads to many neglected microseismic events. Thus detection of weak signal is a vital step in this stage. A typical 1.5s record (8631–8632.5s after the beginning of fracturing) with horizontal components H1 and H2, vertical component V is shown in Fig. 8. As we can see from the initial data, the microseismic record is very noisy and the background noise masks the useful signals. The relative strong S wave is visible in the V component record. However, the events are difficult to follow in both H1 and H2 components. The frequency band of the perforation signal ranges from 0Hz to 500Hz. The frequency band of the observed microseismic signals and background noise ranges from 0Hz to 350Hz and from 0Hz to 900Hz, respectively. Due to the impact of industrial electricity, there are low-frequency interferences in several traces. The S/N of the initial dataset is approximately −14.5942 dB.

3-C microseismic data. (**a**) Horizontal components (H1). (**b**) Horizontal components (H2). (**c**) Vertical component (V).

We then decompose the initial data into five morphological scale components. Fig. 9 shows the results after the proposed orthogonalized morphological reconstruction approach. It can be observed that the events are much more clear than that in the raw data. We can easily follow the coherent energy in all H1 (Fig. 9(a)), H2 (Fig. 9(b)), and V (Fig. 9(c)) component records. The S/N of the denoised result is approximately 4.0927 dB. As a comparison, the traditional morphological reconstruction approach is applied to this example. Similarly, the 2nd–4th scale components are used to reconstruct both H1 and H2 components, and the 2nd–5th scale components are used to reconstruct V components. The weighting coefficients are chosen manually as $\mathrm{(1,1,0.5)}$, $\mathrm{(1,1,0.5)}$ and $\mathrm{(1,1,1,1)}$, respectively.

Orthogonalized morphological reconstruction results of (**a**) horizontal components (H1) by 2nd–4th scales components, (**b**) horizontal components (H2) by 2nd–4th scales components, and (**c**) vertical component (V) by 2nd–5th scales components. **...**

The results using the traditional morphological reconstruction approach are shown in Fig. 10. The events are more visible than the initial data, but the proposed approach performs better. The S/N of the denoised result is approximately −3.1015 dB. In order to avoid manually choosing the weighting coefficient σ_{k}, a varimax norm based morphological reconstruction approach can be used, in which the σ_{k} is defined as:

$${\mathrm{\sigma}}_{k}=\mathrm{1/}nor{m}_{v}({\mathbf{c}}_{k}),$$

31

where *n**o**r**m*_{v}( ⋅ ) is the varimax norm^{50}. The results are shown in Fig. 11. The reconstructions of the H1 (Fig. 11(a)) and V (Fig. 11(c)) components are acceptable, but the reconstruction of H2 (Fig. 11(a)) component is unsatisfied. The event is still hardly detected in the H2 component. The S/N of the denoised result is approximately −5.1291 dB.

Traditional morphological reconstruction results of (**a**) horizontal components (H1) by 2nd–4th scales components, (**b**) horizontal components (H2) by 2nd–4th scales components, and (**c**) vertical component (V) by 2nd–5th scales components. **...**

Varimax norm based morphological reconstruction results of (**a**) horizontal components (H1) by 2nd–4th scales components, (**b**) horizontal components (H2) by 2nd–4th scales components, and (**c**) vertical component (V) by 2nd–5th scales **...**

Event location is an important step in the processing of microseismic data. For further evaluation of denoising performance by our proposed and other competitive approaches, we pick the events in each processed results. An accurate time-picking corresponds to a good weak signal detecting performance. The automatic events detection algorithm is chosen as the well-known STA/LTA filter^{51}. In this test, the energy is used as characteristic function (CF) in STA/LTA filter. In the following figures the time picks are represented with red asterisks. Figure 12 shows the picked arrival times for the noisy 3-C records. Figures 13–15 show the results using the proposed method, median filtering and singular spectrum analysis (SSA) method, respectively. As we can see from this test, because of the strong background noise, the microseismic events are hard to pick. We can observe from Fig. 12 that STA/LTA filter is triggered at incorrect time for many traces. It is obvious that after using the proposed orthogonalized morphological reconstruction approach, the events become much more clear and easier to pick than others, which indicates the superior performance of our proposed approach.

Arrival picking results using median filtering. (**a**) Horizontal component (H1). (**b**) horizontal component (H2). (**c**) Vertical component (V).

Arrival picking of the initial data. (**a**) Horizontal component (H1). (**b**) horizontal component (H2). (**c**) Vertical component (V).

Arrival picking of orthogonalized morphological reconstructions. (**a**) Horizontal component (H1). (**b**) horizontal component (H2). (**c**) Vertical component (V).

Arrival picking results using SSA. (**a**) Horizontal component (H1). (**b**) horizontal component (H2). (**c**) Vertical component (V).

We apply the STA/LTA algorithm to a longer duration of recorded data and denoised data. The experimental results as presented in Table 1 show that more events have been detected after using the proposed denoising approach than using other methods. We use the detected events in the denoised data by the proposed approach to locate the sources. We use Geiger’s approach^{52} to obtain the location. One can find the details of this approach in^{53}. The calculation of travel-time is based on the principle of ray tracing. The results of locating is shown in Fig. 16. The black curve line denotes the trajectory of the fracturing well. The blue circle denotes the position of perforation. The green asterisks denote the locations of the microseismic events.

Comparison of events detection in recorded data and denoised data after using different denoising approaches.

We have proposed a novel denoising method based on mathematical morphological decomposition. We introduce an orthogonalization operator into the process of reconstruction, which can impel the reconstruction of weak signal. We give detailed mathematical introduction of the new method and connect it with several well-known methods and mathematical models. The most striking difference between the proposed and traditional methods is that the core calculations in the proposed method are based on logical operation and set theory. Synthetic and real data examples demonstrate its superior performance compared with the competing alternative approaches. The detected weak signals make the microseismic monitoring feasible in severe environment where the recorded data is extremely noisy and microseismic signals are very weak. The proposed orthogonalized morphological reconstruction method belongs to a class of single-channel techniques and does not require array data. It can be used not only in microseismic monitoring, but also in other type of seismic data (active source or earthquake data), and in other real world applications, e.g., image processing and signal processing, large-scale earthquake data processing and inversion. The proposed method is promising for a wide research community and industrial applications.

The authors would like to thank Yimin Yuan and Sergey Fomel for the constructive discussions. This work is supported by the National Basic Research Program of China (973 Program), grant NO: 2013CB228602. W.H. would like to thank the China Scholarship Council for the financial support.

Author Contributions

W.H. conducted the numerical experiments. R.W. supervised the work. W.H. and H.L. processed the field data. W.H. and Y.C. analyzed the results. W.H., R.W., and Y.C. wrote the manuscript. All authors reviewed the manuscript.

The authors declare that they have no competing interests.

**Electronic supplementary material**

**Supplementary information** accompanies this paper at 10.1038/s41598-017-09711-2.

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

1. Shapiro, S., Dinske, C. & Rothert, E. Hydraulic-fracturing controlled dynamics of microseismic clouds. *Geophysical Research Letters***33** (2006).

2. Grigoli F, et al. Automated microseismic event location using master-event waveform stacking. Scientific Reports. 2016;6:1–13. doi: 10.1038/srep25744. [PMC free article] [PubMed] [Cross Ref]

3. Warpinski NR. Analytic crack solutions for tilt fields around hydraulic fractures. Journal of Geophysical Research: Solid Earth. 2000;105:23463–23478. doi: 10.1029/2000JB900211. [Cross Ref]

4. Parotidis, M., Shapiro, S. & Rothert, E. Back front of seismicity induced after termination of borehole fluid injection. *Geophysical Research Letters***31** (2004).

5. House L. Locating microearthquakes induced by hydraulic fracturing in crystalline rock. Geophysical Research Letters. 1987;14:919–921. doi: 10.1029/GL014i009p00919. [Cross Ref]

6. Phillips W, Fairbanks T, Rutledge J, Anderson D. Induced microearthquake patterns and oil-producing fracture systems in the austin chalk. Tectonophysics. 1998;289:153–169. doi: 10.1016/S0040-1951(97)00313-2. [Cross Ref]

7. Rutledge JT, Phillips WS, Schuessler BK. Reservoir characterization using oil-production-induced microseismicity, clinton county, kentucky. Tectonophysics. 1998;289:129–152. doi: 10.1016/S0040-1951(97)00312-0. [Cross Ref]

8. Rutledge JT, Phillips WS, Mayerhofer M. Faulting induced by forced fluid injection and fluid flow forced by faulting: An interpretation of hydraulic-fracture microseismicity, carthage cotton valley gas field, texas. Bulletin of the Seismological Society of America. 2004;94:1817–1830. doi: 10.1785/012003257. [Cross Ref]

9. Fischer, T., Hainzl, S., Eisner, L., Shapiro, S. & Le Calvez, J. Microseismic signatures of hydraulic fracture growth in sediment formations: Observations and modeling. *Journal of Geophysical Research: Solid Earth***113** (2008).

10. Šlen, J., Hill, D. P., Eisner, L. & Cornet, F. H. Non–double-couple mechanisms of microearthquakes induced by hydraulic fracturing. *Journal of Geophysical Research: Solid Earth***114** (2009).

11. Flewelling SA, Tymchak MP, Warpinski N. Hydraulic fracture height limits and fault interactions in tight oil and gas formations. Geophysical Research Letters. 2013;40:3602–3606. doi: 10.1002/grl.50707. [Cross Ref]

12. Dong L, Li X. A microseismic/acoustic emission source location method using arrival times of ps waves for unknown velocity system. International Journal of Distributed Sensor Networks. 2013;9:307489. doi: 10.1155/2013/307489. [Cross Ref]

13. Dong, L., Li, X. & Xie, G. Nonlinear methodologies for identifying seismic event and nuclear explosion using random forest, support vector machine, and naive bayes classification. In *Abstract and Applied Analysis*, vol. 2014 (Hindawi Publishing Corporation, 2014).

14. Dong, L.-J., Wesseloo, J., Potvin, Y. & Li, X.-B. Discriminant models of blasts and seismic events in mine seismology. *International Journal of Rock Mechanics and Mining Sciences* 282–291 (2016).

15. Dong L, Wesseloo J, Potvin Y, Li X. Discrimination of mine seismic events and blasts using the fisher classifier, naive bayesian classifier and logistic regression. Rock Mechanics and Rock Engineering. 2016;49:183–211. doi: 10.1007/s00603-015-0733-y. [Cross Ref]

16. Maxwell S, Rutledge J, Jones R, Fehler M. Petroleum reservoir characterization using downhole microseismic monitoring. Geophysics. 2010;75:75A129–75A137. doi: 10.1190/1.3477966. [Cross Ref]

17. Bolton H, Masters G. Travel times of p and s from the global digital seismic networks: Implications for the relative variation of p and s velocity in the mantle. Journal of Geophysical Research: Solid Earth. 2001;106:13527–13540. doi: 10.1029/2000JB900378. [Cross Ref]

18. Reyes-Montes, J., Rietbrock, A., Collins, D. & Young, R. Relative location of excavation induced microseismicity at the underground research laboratory (aecl, canada) using surveyed reference events. *Geophysical research letters***32** (2005).

19. Michelet, S. & Toksöz, M. N. Fracture mapping in the soultz-sous-forêts geothermal field using microearthquake locations. *Journal of Geophysical Research: Solid Earth***112** (2007).

20. Chouet, B. *et al*. Source mechanisms of explosions at stromboli volcano, italy, determined from moment-tensor inversions of very-long-period data. *Journal of Geophysical Research: Solid Earth***108** (2003).

21. Li H, Wang R, Cao S, Chen Y, Huang W. A method for low-frequency noise suppression based on mathematical morphology in microseismic monitoring. Geophysics. 2016;81:V159–V167. doi: 10.1190/geo2015-0222.1. [Cross Ref]

22. Chen Y, Ma J. Random noise attenuation by *f-x* empirical mode decomposition predictive filtering. Geophysics. 2014;79:V81–V91. doi: 10.1190/geo2013-0080.1. [Cross Ref]

23. Justusson, B. Median filtering: Statistical properties. In Two-Dimensional Digital Signal Prcessing II, 161–196 (Springer, 1981).

24. Neelamani R, Baumstein AI, Gillard DG, Hadidi MT, Soroka WL. Coherent and random noise attenuation using the curvelet transform. The Leading Edge. 2008;27:240–248. doi: 10.1190/1.2840373. [Cross Ref]

25. Chen Y, Fomel S, Hu J. Iterative deblending of simultaneous-source seismic data using seislet-domain shaping regularization. Geophysics. 2014;79:V179–V189. doi: 10.1190/geo2013-0449.1. [Cross Ref]

26. Mousavi, S. M. & Langston, C. A. Hybrid seismic denoising using higher-order statistics and improved wavelet block thresholding. *Bulletin of the Seismological Society of America* (2016).

27. Chen Y, et al. Simultaneous denoising and reconstruction of 5-D seismic data via damped rank-reduction method. Geophysical Journal International. 2016;206:1695–1717. doi: 10.1093/gji/ggw230. [Cross Ref]

28. Huang W, Wang R, Yuan Y, Gan S, Chen Y. Signal extraction using randomized-order multichannel singular spectrum analysis. Geophysics. 2016;82:V69–V84. doi: 10.1190/geo2015-0708.1. [Cross Ref]

29. Kong S, Phinney RA, Roy-Chowdhury K. A nonlinear signal detector for enhancement of noisy seismic record sections. Geophysics. 1985;50:539–550. doi: 10.1190/1.1441931. [Cross Ref]

30. Schimmel M, Paulssen H. Noise reduction and detection of weak, coherent signals through phase-weighted stacks. Geophysical Journal International. 1997;130:497–505. doi: 10.1111/j.1365-246X.1997.tb05664.x. [Cross Ref]

31. Gibbons SJ, Ringdal F. The detection of low magnitude seismic events using array-based waveform correlation. Geophysical Journal International. 2006;165:149–166. doi: 10.1111/j.1365-246X.2006.02865.x. [Cross Ref]

32. Mousavi SM, Langston CA, Horton SP. Automatic microseismic denoising and onset detection using the synchrosqueezed continuous wavelet transform. Geophysics. 2016;81:V341–V355. doi: 10.1190/geo2015-0598.1. [Cross Ref]

33. Mousavi SM, Langston CA. Adaptive noise estimation and suppression for improving microseismic event detection. Journal of Applied Geophysics. 2016;132:116–124. doi: 10.1016/j.jappgeo.2016.06.008. [Cross Ref]

34. Matheron, G. *Random sets and integral geometry* (John Wiley & Sons, 1975).

35. Serra, J. Image analysis and mathematical morphology, v. 1 (Academic press, 1982).

36. Koskinen, L., Astola, J. T. & Neuvo, Y. A. Soft morphological filters. In *San Diego,’91, San Diego, CA*, 262–270 (International Society for Optics and Photonics, 1991).

37. Sinha D, Dougherty ER. Fuzzy mathematical morphology. Journal of Visual Communication and Image Representation. 1992;3:286–302. doi: 10.1016/1047-3203(92)90024-N. [Cross Ref]

38. Wang R, Li Q, Zhang M. Application of multi-scaled morphology in denoising seismic data. Applied Geophysics. 2008;5:197–203. doi: 10.1007/s11770-008-0033-3. [Cross Ref]

39. Fomel S. Shaping regularization in geophysical-estimation problems. Geophysics. 2007;72:R29–R36. doi: 10.1190/1.2433716. [Cross Ref]

40. Tikhonov A. Solution of incorrectly formulated problems and the regularization method. Soviet Math. Dokl. 1963;5:1035–1038.

41. Maragos P. Morphological systems: slope transforms and max-min difference and differential equations. Signal Processing. 1994;38:57–77. doi: 10.1016/0165-1684(94)90057-4. [Cross Ref]

42. Van Huffel S. Enhanced resolution based on minimum variance estimation and exponential data modeling. Signal processing. 1993;33:333–355. doi: 10.1016/0165-1684(93)90130-3. [Cross Ref]

43. Chen Y, Fomel S. Random noise attenuation using local signal-and-noise orthogonalization. Geophysics. 2015;80:WD1–WD9. doi: 10.1190/geo2014-0227.1. [Cross Ref]

44. Huang W, Wang R, Chen Y, Li H, Gan S. Damped multichannel singular spectrum analysis for 3D random noise attenuation. Geophysics. 2016;81:V261–V270. doi: 10.1190/geo2015-0264.1. [Cross Ref]

45. Huang W, Wang R, Chen X, Chen Y. Double least-squares projections method for signal estimation. IEEE Transactions on Geoscience & Remote Sensing. 2017;55:4111–4129. doi: 10.1109/TGRS.2017.2688420. [Cross Ref]

46. Fomel S. Adaptive multiple subtraction using regularized nonstationary regression. Geophysics. 2008;74:V25–V33. doi: 10.1190/1.3043447. [Cross Ref]

47. Hestenes, M. R. & Stiefel, E. *Methods of conjugate gradients for solving linear systems*, vol. 49 (NBS, 1952).

48. Chen, Y. *et al*. Empirical low-rank approximation for seismic noise attenuation. *IEEE Transactions on Geoscience and Remote Sensing* (2017).

49. Stockwell RG, Mansinha L, Lowe R. Localization of the complex spectrum: the s transform. IEEE transactions on signal processing. 1996;44:998–1001. doi: 10.1109/78.492555. [Cross Ref]

50. Wiggins RA. Minimum entropy deconvolution. Geoexploration. 1978;16:21–35. doi: 10.1016/0016-7142(78)90005-4. [Cross Ref]

51. Allen RV. Automatic earthquake recognition and timing from single traces. Bulletin of the Seismological Society of America. 1978;68:1521–1532.

52. Geiger L. Probability method for the determination of earthquake epicenters from the arrival time only. Bull. St. Louis Univ. 1912;8:56–71.

53. Ge M. Analysis of source location algorithms part ii: iterative methods. Journal of Acoustic Emission. 2003;21:29–51.

Articles from Scientific Reports are provided here courtesy of **Nature Publishing Group**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |