TY - JOUR
AU - Chen, Yangkang
AB - Summary Microseismic signal is typically weak compared with the strong background noise. In order to effectively detect the weak signal in microseismic data, we propose a mathematical morphology based approach. We decompose the initial data into several morphological multiscale components. For detection of weak signal, a non-stationary weighting operator is proposed and introduced into the process of reconstruction of data by morphological multiscale components. The non-stationary weighting operator can be obtained by solving an inversion problem. The regularized non-stationary method can be understood as a non-stationary matching filtering method, where the matching filter has the same size as the data to be filtered. In this paper, we provide detailed algorithmic descriptions and analysis. The detailed algorithm framework, parameter selection and computational issue for the regularized non-stationary morphological reconstruction (RNMR) method are presented. We validate the presented method through a comprehensive analysis through different data examples. We first test the proposed technique using a synthetic data set. Then the proposed technique is applied to a field project, where the signals induced from hydraulic fracturing are recorded by 12 three-component geophones in a monitoring well. The result demonstrates that the RNMR can improve the detectability of the weak microseismic signals. Using the processed data, the short-term-average over long-term average picking algorithm and Geiger’s method are applied to obtain new locations of microseismic events. In addition, we show that the proposed RNMR method can be used not only in microseismic data but also in reflection seismic data to detect the weak signal. We also discussed the extension of RNMR from 1-D to 2-D or a higher dimensional version. Image processing, Earthquake source observation, Seismic Noise INTRODUCTION It has already been recognized that microseismic monitoring has a significant potential to characterize physical processes related to fluid injections and extractions in hydrocarbon and geothermal reservoirs (Shapiro et al. 2006; Xia et al. 2013; Zhang et al. 2015; Mousavi et al. 2016a,b; Cheng et al. 2016; Ebrahimi et al. 2016; Chen et al. 2017; Guo & Mcmechan 2017; Guo et al. 2017; Lu & Feng 2017; Mousavi & Langston 2017; Shahin et al. 2017; Tian et al. 2017; Zhang et al. 2017). Localization of the associated microseismic events enables mapping and calibrating stimulations. This technique has been widely studied and applied in petroleum and gas exploration (House 1987; Phillips et al. 1998; Rutledge et al. 1998, 2004; Gibbons & Ringdal 2006a; Shapiro et al. 2006; Gibbons et al. 2008; Fischer et al. 2008; Šílený et al. 2009; Flewelling et al. 2013; Chen 2018). Specifically, Pearson (1981) studied the relationship between microseismicity and high pore pressures during hydraulic stimulation experiments, using a 1-D diffusion model. Fehler et al. (1998) used a parametrization scheme in which the velocities of each block are allowed to change from the background velocity only after a threshold number of microearthquakes have occurred in the block, in order to improve image of the velocity structure. Reshetnikov et al. (2010) developed a two-step passive seismic imaging approach, in which the hypocentre of the microseismic event is first located, and then the reflections within the recorded wavefield are processed by using a directional migration algorithm. Hansen & Schmandt (2015) combined a reverse-time imaging method and short-term-average over long-term average technique to automatically detect and locate microseismicity. Hajati et al. (2015) analysed the interevent time distribution of hydraulic-fracturing-induced seismicity and utilize a universal statistical process to describe the distribution of hydraulic-fracturing-induced events in time. 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 (Maxwell et al. 2010). In this paper, we use the term ‘weak signal’ to refer to the signal with a very low signal-to-noise ratio (SNR). The weak signal is hardly to be observed directly from the profile. A poor SNR can lead to unauthentic arrival time-picks (Bolton & Masters 2001) and localization of microseismic events (Reyes-Montes et al. 2005). Moreover, weak signal is easily masked, resulting in loss of microseismic events. All of these will negatively affect the performance when imaging those microseismic activities and fractures (Michelet & Toksöz 2007) and inverting source mechanisms (Chouet et al. 2003) for better fracture characterization. Therefore the detection of weak signal and improvement of the SNR are necessary prerequisites in microseismic monitoring. Over the past few decades, many techniques coming from various fields have been proposed to handle the noise issues, including median filtering (MR; Justusson 1981; Bai & Wu 2017); mean filtering (Xie et al. 2015b, 2016); decomposition and reconstruction (Han & Mirko 2015; Chen et al. 2016c, 2016b; Liu et al. 2016); various kinds of mathematical transform based approaches (Huang et al. 2017c; Mortazavi et al. 2017), for example, curvelet transform (Neelamani et al. 2008; Wang & Wang 2017), dreamlet transform (Wang et al. 2015a), synchrosqueezing wavelet transform (Xie et al. 2015a) and seislet transform (Chen & Fomel 2015b; Gan et al. 2015a,b; Chen 2016); matrix completion based approach (Chen et al. 2016e, 2016d; Zhang et al. 2016c, 2016b, 2016a; Huang et al. 2016a,b); various signal detectors, for example, spatial coherence and slowness (Kong et al. 1985), instantaneous phase (Schimmel & Paulssen 1997) and array-based waveform correlation (Gibbons & Ringdal 2006b) based signal detectors. In microseismic monitoring, the most commonly used method for reducing background noise and detecting weak signal is frequency filtering (Li et al. 2016b). The frequency filtering utilizes the frequency differences between signal and noise to separate them. However, its application is limited when the signal and noise overlap in the frequency domain. Wang (2005) first introduced mathematical morphological filtering (MMF) in seismic data processing. In the work presented by Wang (2005), the morphological difference between signal and noise was first proposed in seismics, and used to remove abnormal amplitudes. MMF tends to simplify image data preserving their essential shape characteristics and eliminating irrelevancies. Unlike traditional methods in seismic data processing, the basis of MMF are logical operation and set theory, which can provide us a method to process signal over the complete frequency bandwidth, improving the SNR and maintaining the resolution (Wang et al. 2008). Mathematical morphology (Matheron 1975; Serra 1982; Wang et al. 2008; Huang et al. 2017b) 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. Matheron (1975) and Serra (1982) set up the foundation of mathematical morphology. Then, a number of improvements and modifications for this method were proposed such as the soft mathematical morphology (Koskinen et al. 1991) and the generalized mathematical morphology (Sinha & Dougherty 1992). Later, this method was rapidly developed and widely applied in seismic data processing. For example, Li et al. (2016a) 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. Huang et al. (2017d) developed MMF from the time direction to the spatial direction, and used it to suppress coherent noise, with the combination of a local slope scanning algorithm. Zhou et al. (2016) combined MMF and a mesh scanning algorithm, and utilized it to attenuate the diffraction noise in marine surveys. The multiscale morphological decomposition can decompose object signal into different components with different morphological scales, by repeated MMFs with a series of structuring elements (SE), which have the same shape but different scales. Wang et al. (2008) first used morphological decomposition and reconstruction to suppress ground roll in land seismic data, by abandoning those components with greater morphological scales in the process of reconstruction, in which, most energy of the ground roll is included. Chen et al. (2015) applied the multiscale morphology theory into the detection of the weak signal in microseismic monitoring. Yuan et al. (2016) proposed a Varimax norm (Wiggins 1978) based approach to adaptively choose weight factors of reconstruction in multiscale morphological decomposition, and applied it to detect weak signal in thin interbedded reservoirs. In this paper, a new regularized non-stationary morphological reconstruction (RNMR) algorithm is proposed, as a pre-processing technique, for the detection of microseismic weak signal. In the new RNMR approach, a non-stationary weighting operator is introduced as a projection operator that projects the input signal on a subspace spanned by several selected morphological basis vectors. It can be also understood as a non-stationary counterpart of the well-known matching filter, where the length of the matching filter is the same as the input signal. Unlike the former morphological reconstruction (MR) approaches, which are manual and subjective (Wang et al. 2008; Chen et al. 2015), or empirical, as used in Yuan et al. (2016), the new RNMR method solves the problem of reconstruction of signal using an inversion framework. A regularization, or in other words, penalty term is necessary to constrain the solution space and guarantee the stability of the inversion problem. One of the most classical and commonly used regularization approaches is Tikhonov’s regularization (Tikhonov 1963). Although effective, as Fomel (2007b) indicated, it is difficult to find and specify an appropriate regularization operator. We choose the more convenient shaping regularization approach (Fomel 2007b) for solving the inversion problem in the RNMR algorithm. In fact, the shaping regularization approach has been widely used in geophysical problems, such as local seismic attributes (Fomel 2007a), time-frequency analysis (Liu et al. 2011), non-stationary decomposition (Wu et al. 2016), noise attenuation (Huang et al. 2017a), interpolation (Zu et al. 2016b), and deblending (Chen et al. 2014; Zu et al. 2015, 2016a, 2017). A complete and detailed algorithm workflow for RNMR will be presented. We first use the synthetic data experiments to test the performance of the proposed RNMR approach. Then the proposed approach is applied to a real three-component (3-C) microseismic data set. Compared with the former MR approaches, bandpass filtering, MR (Justusson 1981), and singular spectrum analysis (SSA; Trickett 2008; Huang et al. 2016a) approach, the proposed approach demonstrates a superior performance. THEORY Theory of morphology Let d be a seismic trace and $${\bf f}\subseteq \mathbb {R}$$ be a set of amplitude values. The value of a point t in d is represented by d(t) ∈ f. The morphological dilation ϕb(d) and erosion φb(d) are the morphological operations that transform d with the SE b(τ) as (Maragos 1994):   \begin{equation} \phi _b({\bf d})=\bigvee _\tau d(t-\tau )+b(\tau ), \end{equation} (1)  \begin{equation} \varphi _b({\bf d})=\bigwedge _\tau d(t+\tau )-b(\tau ), \end{equation} (2)where $$\bigvee$$ denotes supremum, and $$\bigwedge$$ denotes infimum. The SE b(τ) is a small section of a symmetric and convex signal (function) with a specific shape. The detailed discussion of SE will be given in the following section. 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 (dilation) and morphological dilation (erosion) forms the morphological opening χb(d) (closing ψb(d)) as:   \begin{equation} \chi _b({\bf d})=\phi _b(\varphi _b({\bf d})), \end{equation} (3)  \begin{equation} \psi _b({\bf d})=\varphi _b(\phi _b({\bf d})). \end{equation} (4) Consider $$\lbrace \chi _{b_k}({\bf d})\rbrace$$, k ∈ [1, K] and $$\lbrace \psi _{b_k}({\bf d})\rbrace$$, k ∈ [1, K], two indexed families of morphological opening and closing, respectively. Typically, the index k denotes the morphological scale. To avoid notational clutter, let $$f_{b_k}=(\chi _{b_k}\psi _{b_k}+\psi _{b_k}\chi _{b_k})/2$$. d can be represented by an additive decomposition with K+1 scales as follow:   \begin{eqnarray} {\bf d}&=&f_{b_0}({\bf d})-f_{b_1}f_{b_0}({\bf d})+f_{b_1}f_{b_0}({\bf d}) \nonumber\\ &=&f_{b_0}({\bf d})-\!f_{b_1}f_{b_0}({\bf d})+f_{b_1}f_{b_0}({\bf d})-\!f_{b_2}f_{b_1}f_{b_0}({\bf d})+f_{b_2}f_{b_1}f_{b_0}({\bf d})\nonumber\\ &&......\nonumber\\ &=&\sum _{k=1}^K[f_{b_{k-1}}...f_{b_{1}}f_{b_{0}}({\bf d})-f_{b_{k}}...f_{b_{1}}f_{b_{0}}({\bf d})]+f_{b_{k}}...f_{b_{1}}f_{b_{0}}({\bf d}), \end{eqnarray} (5) where $$f_{b_0}({\bf d})={\bf d}$$. Eq. (5) refers to multiscale morphological decomposition (Wang et al. 2008). $$f_{b_k}$$ actually refers to MMF (Wang 2005; Zhou et al. 2016). The multiscale morphological decomposition is essentially the repetitive MMF with a series of SE. For convenience, let:   \begin{eqnarray} {\bf c}_k=\left\lbrace \begin{array}{ll}f_{b_{k-1}}...f_{b_{1}}f_{b_{0}}({\bf d})-f_{b_{k}}...f_{b_{1}}f_{b_{0}}({\bf d}), &\quad k\in [1,K]\hbox{,} \\ f_{b_{k}}...f_{b_{1}}f_{b_{0}}({\bf d}), &\quad k=K+1\hbox{,} \end{array} \right. \end{eqnarray} (6)where ck, k ∈ [1, K + 1] are called the morphological multiscale components (Wang et al. 2008). The value of a point t in ck is represented by ck(t) ∈ f. Our experience is that 3-10 decomposed components are appropriate for most seismic data sets, taking the compromise between efficiency and effectiveness into consideration. Algorithm 1 gives a detailed description of the multiscale morphological decomposition. The nature of the multiscale morphological decomposition is to decompose the discrete data set d into a series of primary subsets ck (Wang et al. 2008), which satisfies that:   \begin{equation} {\bf d}=\bigcup _{k=1}^{K+1}{\bf c}_{k}, \end{equation} (7)  \begin{equation} \emptyset =\bigcap _{k=1}^{K+1}{\bf c}_{k}, \end{equation} (8)where ∅ denotes empty set. The reconstruction of data by ck can be typically represented as (Wang et al. 2008):   \begin{equation} \sum _{{\bf c}_k\in {\bf E}}[w_k{\bf c}_k]=\tilde{{\bf d}}, \end{equation} (9)where E is a subset of {ck}k ∈ [1, K], that E⊆{ck}k ∈ [1, K]. Constant wk ∈ [0, 1] is the weighting coefficient which can control energy from difference scale components. This decomposition allows for full reconstruction of the original data, when E = {ck}k ∈ [1, K] and ωk ≡ 1. Algorithm 1 View largeDownload slide Multiscale morphological decomposition Algorithm 1 View largeDownload slide Multiscale morphological decomposition Structuring element The SE is key in the morphological decomposition and reconstruction, which 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 straight line, a semicircle or a triangle, defined respectively as:   \begin{equation} b_{st}(\tau )=A, \end{equation} (10)  \begin{equation} b_{se}(\tau )=A \sqrt{1-(\tau /W)^2}, \end{equation} (11)  \begin{equation} b_{tr}(\tau )=A (1-|\tau /W|), \end{equation} (12)where τ ∈ [ − W, W], W > 0. bst, bse and btr represent the straight line, semicircle and triangle SEs, respectively. A and W denote the height and width. Fig. 1 shows a demonstration of the three types of SE. The SE with different combination of the three parameters has different morphological scale. When the shape of an SE is fixed, its scale increases when the height decreases, or the width increases. On the contrary, its scale decreases when the height increases, or the width decreases. The comparison of scale among the three shapes is as follows:   \begin{eqnarray} Scale(straightline)>Scale(semicircle)>Scale(triangle). \nonumber\\ \end{eqnarray} (13) Figure 1. View largeDownload slide Three types of SE. (a) Straight line SE. (b) Semicircle SE. (c) Triangle SE. Figure 1. View largeDownload slide Three types of SE. (a) Straight line SE. (b) Semicircle SE. (c) Triangle SE. Repetitive MMF with a series of SE with different scales (i.e. multiscale morphological decomposition) can decompose the input data into different components and obtain the different morphological information of the input data. Generally, for a specific morphological decomposition, the shapes of the SE family bK are fixed but their heights and widths increase gradually. The increase rate determines the fineness of decomposition. A commonly used strategy to create the SE family bK is using self-morphological dilation. Given a seed SE b1, the ith SE bi can be produced by i − 1 times self-morphological dilation written as:   \begin{equation} {\bf b}_i=\underbrace{\phi _{b_1}(...\phi _{b_1}(\phi _{b_1}}_{i-1\ {\rm times}}({\bf b}_1))). \end{equation} (14) Fig. 2 shows a demonstration of applying the self-morphological dilation of the three types of SE four times. The smallest SEs (b1) are the seed SEs. The seed SEs are self-dilated four times and produce b2, b3, b4, and b5. Fig. 3 gives an example of a seven scale morphological decomposition of a Ricker wavelet with various scales of semicircle SE. The 1st trace (scale 0) is the initial wavelet. The second to eighth traces are the seven scale components. The choice of the shape of SE is related to the input signal and the processing target. For example, the straight line SE can be used in seismic coherent noise suppression (Huang et al. 2017d), where MMF is implemented by sliding the straight line SE along the trajectory of coherent noise (spatial direction). The straight line SE has a better performance in separating the energy of the coherent noise and others in the direction along the trajectory of coherent noise. The semicircle SE is generally used in the temporal direction of each single trace to detect the morphological information of signal (Wang 2005; Li et al. 2016a) because its shape can better match seismic waveform and obtain smoother results compared with the straight line SE. Therefore the SE used in this research is semicircular. The triangle SE is rarely used in seismic data processing so far, but it is widely used in image processing such as corner detection (Sobania & Evans 2005). Figure 2. View largeDownload slide Demonstration of the self-morphological dilation of the three types of SE. (a) Self-morphological dilation of straight line SE. (b) Self-morphological dilation of semicircle SE. (c) Self-morphological dilation of triangle SE. The smallest SEs (b1) in each panel are the seed SEs. Figure 2. View largeDownload slide Demonstration of the self-morphological dilation of the three types of SE. (a) Self-morphological dilation of straight line SE. (b) Self-morphological dilation of semicircle SE. (c) Self-morphological dilation of triangle SE. The smallest SEs (b1) in each panel are the seed SEs. Figure 3. View largeDownload slide Multiscale morphological decomposition of a Ricker wavelet. The first trace (scale 0) is the initial wavelet. The second to eighth traces are the seven scale components. Figure 3. View largeDownload slide Multiscale morphological decomposition of a Ricker wavelet. The first trace (scale 0) is the initial wavelet. The second to eighth traces are the seven scale components. Regularized non-stationary morphological reconstruction (RNMR) In traditional MR (Wang et al. 2008; Chen et al. 2015), the weighting coefficient wk is chosen manually, which makes the reconstruction subjective. In addition, it is difficult to choose an appropriate weighting coefficient, and the process of choosing costs a lot of time and human effort. For weak signal detection, we define a new MR as an inversion problem:   \begin{equation} \arg \min _{w_k(t)}\sum _{{\bf c}_k\in {\bf E}}[\parallel w_k(t)c_k(t)-d(t)\parallel _2^2+\mathscr {R}(w_k(t))], \end{equation} (15)where ∥ · ∥2 represents the L2 norm of a vector, and $$\mathscr {R}$$ represents the regularization operator. Observe that eq. (15) transfers the constraint of the weighting coefficient to the regularization term. On the one hand, it loosens the constraint of wk and allow it to change with t, in other words, wk(t) now is a non-stationary weighting operator. On the other hand, it puts the constraint to wk by regularizing it, which maps the model (wk) to the space of admissible functions. Hence, the new MR (eq. 15) is named as RNMR. We rewrite eq. (15) using matrix form:   \begin{equation} \arg \min _{\boldsymbol {W}_k}\sum _{{\bf c}_k\in {\bf E}}[\parallel \boldsymbol {W}_k{\bf c}_k-{\bf d}\parallel _2^2+\mathscr {R}(\boldsymbol {W}_k)], \end{equation} (16)where Wk is a diagonal matrix composed by wk(t): Wk = diag(wk(t)). The RNMR problem (eq. 16) is solved by shaping regularization (Fomel 2007b), which will be discussed in the following section. The estimated signal $$\bar{{\bf d}}$$ by RNMR holds as:   \begin{equation} \sum _{{\bf c}_k\in {\bf E}}[\boldsymbol {W}_k{\bf c}_k]=\bar{{\bf d}}. \end{equation} (17) The proposed non-stationary weighting operator is actually a projection operator that projects the initial data d (a higher-dimensional vector) on the space spanned by several selected morphological basis vectors (a lower-dimensional space). Regularizing non-stationary weighting via shaping regularization Using Wkck = CkWk = Wk · ck, we rewrite eq. (16) as:   \begin{equation} \arg \min _{\boldsymbol {w}_k}\sum _{{\bf c}_k\in {\bf E}}[\parallel {\bf C}_k\boldsymbol {w}_k-{\bf d}\parallel _2^2+\mathscr {R}(\boldsymbol {w}_k)], \end{equation} (18)where Ck is a diagonal matrix composed by ck(t): Ck = diag(ck(t)). Wk is a column vector composed by wk(t): Wk = [wk(t)]T. Now eq. (18) is a regularized linear least-squares problem where Ck is the forward operator, Wk is the model and d is the observed data. It is actually to find a low-dimensional approximation of d, but the model (Wk) is under certain constraints. Tikhonov strategy (Tikhonov 1963) amounts to a particular choice of the regularization operator $$\mathscr {R}(\boldsymbol {w}_k)$$, in which one additionally attempts to minimize the norm of TWk, where T is the regularization operator. Thus eq. (18) can be written as:   \begin{equation} \arg \min _{\boldsymbol {w}_k}\sum _{{\bf C}_k\in {\bf E}}[\parallel {\bf C}_k\boldsymbol {w}_k-{\bf d}\parallel _2^2+\varepsilon ^2\parallel {\bf T}\boldsymbol {w}_k\parallel _2^2], \end{equation} (19)where ε is a scalar scaling parameter. Set the derivatives of the above-mentioned object function with respect to Wk and obtain the solution as:   \begin{equation} \bar{\boldsymbol {w}_k}=({\bf C}_k^T{\bf C}_k+\varepsilon ^2{\bf T}^T{\bf T})^{-1}{\bf C}_k^T{\bf d}. \end{equation} (20)where $$\bar{\boldsymbol {w}_k}$$ is a least-squares estimation of Wk. 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’s perspective (Fomel 2008). Fomel (2007b) proposed using a shaping operator H to constrain the model (Wk), which is easier to choose and the inversion with the shaping regularization is both fast-converging and producing reasonable results at the very first iterations. The shaping operator is controlled by the smoothing radius (i.e. the width of the triangle smoother) (Fomel 2007b). The basic idea of shaping regularization is recognizing smoothing as a fundamental operation. The task of smoothing is to find a model that fits the observed data, but is in a certain sense smoother (Fomel 2007b). In this case, the forward operator is simply the identity operator. Thus the shaping operator H can be represented as (Fomel 2007b):   \begin{equation} \boldsymbol {H}=({\bf I}+\varepsilon ^2{\bf T}^T{\bf T})^{-1}. \end{equation} (21)So,   \begin{equation} \varepsilon ^2{\bf T}^T{\bf T}=\boldsymbol {H}^{-1}-{\bf I}. \end{equation} (22)Combining eqs (20) and (22), we hold the solution of RNMR problem (eq. 18) via shaping regularization as:   \begin{equation} \bar{\boldsymbol {w}_k}=[\lambda ^2{\bf I}+\boldsymbol {H}({\bf C}_k^T{\bf C}_k-\lambda ^2{\bf I})]^{-1}\boldsymbol {H}{\bf C}_k^T{\bf d}. \end{equation} (23)where λ is an introduced parameter controlling the physical dimensionality and enabling fast convergence when inversion is implemented iteratively (Chen et al. 2016e). The shaping operator H can be chosen as certain smoothing operator such as low-pass filtering and Gaussian filtering. In this research, we choose H as a symmetric generalized triangle shaper (Claerbout 1992). The size of the linear system is related to the samples of the input trace. If we assume the input signal d is an n × 1 vector. Ck is an n × n diagonal matrix. H is a 2r − 1 diagonal matrix of size n × n, where r is the smoothing radius. I is a unit matrix of size n × n. The inversion problem in eq. (23) is solved by the conjugate gradient (CG) method (Hestenes & Stiefel 1952), where the shaping operator is split into two matrices to meet the requirement of the CG method. So eq. (23) is rewritten as:   \begin{equation} \bar{\boldsymbol {w}_k}={\bf L}[\lambda ^2{\bf I}+{\bf L}^T({\bf C}_k^T{\bf C}_k-\lambda ^2{\bf I}){\bf L}]^{-1}{\bf L}^T{\bf C}_k^T{\bf d}, \end{equation} (24)where H = LLT. The estimated weak signal s by the RNMR approach can be represented as:   \begin{equation} {\bf s}\approx \sum _{{\bf c}_k\in {\bf E}}[{\bf C}_k\bar{\boldsymbol {w}_k}]. \end{equation} (25)Algorithm 2 gives a detailed description of RNMR. Algorithm 2 View largeDownload slide Regularized non-stationary morphological reconstruction (RNMR) Algorithm 2 View largeDownload slide Regularized non-stationary morphological reconstruction (RNMR) The non-stationary weighting operator in the presented algorithm can be considered as a non-stationary counterpart of the well-known matching filter, where the length of the matching filter is the same as the data (to be filtered). Solving the non-stationary matching filter is extremely challenging because of inverting a highly singular diagonal matrix. Inverting a highly singular diagonal matrix is not computationally difficult, but inverting such a matrix with a reasonable solution is not trivial. A traditional damped least-squares solver would introduce a lot of unphysical artefacts to the solution. The shaping regularization algorithm plays an indispensable role here because it applies a smoothness constraint to the weighting operator, that is, the non-stationary matching filter needs to be locally smooth everywhere, which is reasonable and easy to control via the smoothing radius. The shaping regularization imposes additional constraints to the solution space and guarantees the stability of the inversion problem (Fomel 2007b). The shaping operator provides an explicit mapping of the solution to the space of acceptable solutions. It is embedded in an iterative optimization scheme (the CG algorithm) and allows for a better control on the estimation result. This method can also be viewed as a local orthogonalization algorithm (Chen & Fomel 2015a), where we aim to orthogonalize each scale component. EXAMPLES Synthetic examples We first test the performance of the proposed technique using two synthetic examples. The first example is shown in Fig. 4. As the Ricker wavelet is widely used as the synthetic model in seismological community (Liu et al. 2011; Yang et al. 2014; Jia et al. 2015; Wang et al. 2015a,b; Chen et al. 2016a; Li et al. 2016a; Qu et al. 2016; Chen 2017), we also use it as our synthetic model. The synthetic signal (trace 1) is a Ricker wavelet with 100 Hz dominant frequency and π/2 initial phase. The synthetic noise (trace 2) is broadband Gaussian noise. The first trace is added to the second trace as the input data (trace 3). We can see that the trace 3 is very noisy and the signal is masked by the strong background noise. The seven components using various scales of semicircle SE are shown in traces 4–10. The results of using four bandpass filters (1–10–190–200 Hz, 25–35–165–175 Hz, 50–60–140–150 Hz, 75–85–115–125 Hz) are shown in traces 11–14. The results of using traditional MR, varimax norm based MR (Yuan et al. 2016) and RNMR are presented in traces 15–17. Traces 18–24 demonstrate the corresponding removed noise (differences between trace 3 and traces 11–17). Traces 25–31 demonstrate the corresponding errors (differences between trace 1 and traces 11–17). The varimax norm based MR is proposed to avoid manually choosing the weighting coefficient in traditional MR, where the weighting coefficient wk is chosen automatically as:   \begin{equation} w_k=(\sum _t (c_k(t))^2)^2/\sum _t (c_k(t))^4, \end{equation} (26)where ck(t) is the kth morphological multiscale components. We observe that most energy of the signal is decomposed in the third–seventh scale components (traces 6–10) and trace 10 has a lower SNR compared with traces 6–9. Therefore we choose the reconstruction factor as (1, 1, 1, 1, 0.5) for the traditional MR approach. For the varimax norm based MR and proposed RNMR methods, we use traces 6–10 to reconstruct the signal. The weighting factors in both varimax norm based MR and proposed RNMR are obtained adaptively. It can be seen that although we adjust the frequency band to make the bandpass filtering more and more rigorous, it cannot separate the signal and noise very well because they share the same frequency band. Both traditional and varimax norm based MRs suppress some noise. However, the proposed RNMR approach obtains the best result and reconstructs the signal that is closest to the true signal. For a more detailed comparison, Fig. 5 presents time–frequency spectra of clean data, noisy data and seven processed results, which are calculated by the standard Stockwell transform. It can be observed that the energy of the synthetic signal is concentrated in the area of 0.12–0.17 ms and 0–200 Hz in the time–frequency spectrum (Fig. 5a). Contaminated by Gaussian noise, the time-frequency spectrum (Fig. 5b) becomes noisy. Figs 5(c)–(i) present the results after using the four bandpass filters, the MR, varimax norm based MR and RNMR approaches, respectively. We can observe that the bandpass filter truncates the frequency band of the input data and removes the low and high frequency components. However, it cannot handle the situations where the frequencies of the signal and noise overlap. The proposed RNMR achieves an excellent reconstruction of the signal. We can see that the reconstructed signal (Fig. 5i) is very similar to the initial signal (Fig. 5a). Both conventional and varimax norm based MRs improve the qualities of input data, but their results are still unsatisfactory. Table 1 shows the quantitative analysis, where SNR is used as a statistical indicator to evaluate the performance. The definition of SNR is shown below:   \begin{equation} {\rm SNR}=10\log \frac{\Vert {\bf s}\Vert _2^2}{\Vert {\bf s-d}\Vert _2^2}, \end{equation} (27)where s is the true signal. d denotes the noisy or denoised data. Figure 4. View largeDownload slide The first synthetic example. From left to right, traces 1–2: signal and Gaussian noise; trace 3: signal+Gaussian noise (input data); traces 4–10: seven scale components; traces 11–17: results using four bandpass filters (1–10–190–200 Hz, 25–35–165–175 Hz, 50–60–140–150 Hz, 75–85–115–125 Hz), MR, varimax norm based MR, and RNMR, respectively; traces 18–24: the corresponding removed noise (differences between trace 3 and traces 11–17); traces 25–31: the corresponding errors (differences between trace 1 and traces 11–17). Figure 4. View largeDownload slide The first synthetic example. From left to right, traces 1–2: signal and Gaussian noise; trace 3: signal+Gaussian noise (input data); traces 4–10: seven scale components; traces 11–17: results using four bandpass filters (1–10–190–200 Hz, 25–35–165–175 Hz, 50–60–140–150 Hz, 75–85–115–125 Hz), MR, varimax norm based MR, and RNMR, respectively; traces 18–24: the corresponding removed noise (differences between trace 3 and traces 11–17); traces 25–31: the corresponding errors (differences between trace 1 and traces 11–17). Figure 5. View largeDownload slide Comparison of time–frequency spectra of (a and b) clean and noisy data, (c–f) results after using 1–10–190–200 Hz, 25–35–165–175 Hz, 50–60–140–150 Hz, and 75–85–115–125 Hz bandpass filters, respectively, (g–i) results after using MR, varimax norm based MR, and RNMR, respectively. Figure 5. View largeDownload slide Comparison of time–frequency spectra of (a and b) clean and noisy data, (c–f) results after using 1–10–190–200 Hz, 25–35–165–175 Hz, 50–60–140–150 Hz, and 75–85–115–125 Hz bandpass filters, respectively, (g–i) results after using MR, varimax norm based MR, and RNMR, respectively. Table 1. Quantitative comparison of different approaches for the first synthetic example. BP represents bandpass filtering and Varimax MR represents varimax norm based MR.   Input  BP (1–10–190–200 Hz)  BP (25–35–165–175 Hz)  BP (50–60–140–150 Hz)  SNR  −12.0489  −0.5059  1.2905  2.0579    BP (75–85–115–125 Hz)  MR  Varimax MR  RNMR  SNR  1.4637  −0.7415  −2.5127  11.7489    Input  BP (1–10–190–200 Hz)  BP (25–35–165–175 Hz)  BP (50–60–140–150 Hz)  SNR  −12.0489  −0.5059  1.2905  2.0579    BP (75–85–115–125 Hz)  MR  Varimax MR  RNMR  SNR  1.4637  −0.7415  −2.5127  11.7489  View Large Fig. 6 presents the second synthetic example. The synthetic signal (trace 1) is the same to that in the first synthetic example. The added background noise consists of Gaussian noise (trace 2) and limited-band random noise (trace 3). Traces 2 and 3 are added to trace 1 as the input data (trace 4). Similarly, traces 5–11 show the seven scale morphological components using various scales of semicircle SE, traces 12–15 show the bandpass filtered results (1–10–190–200 Hz, 25–35–165–175 Hz, 50–60–140–150 Hz, 75–85–115–125 Hz), traces 16–18 show the three different morphology based reconstructions, and traces 19–25 and 26–32 show the corresponding removed noise and errors, respectively. In this experiment, traces 7–10 are used for the three different morphology based reconstructions, taking the compromise of preserving signal and suppressing noise into consideration. For the conventional MR, we choose weighting factor as (1, 1, 1, 1). The added limited-band random noise poses some trouble for the detection of signal. The performances of all the methods are worse than those in the first example because of the existence of the limited-band noise, however, the proposed RNMR shows a superior robustness compared with others. The reconstructed signal by the RNMR approach (trace 18) is more detectable than before (trace 4). Fig. 7 presents the time-frequency spectra of the clean, noisy and processed data. It is clear that the strong limited-band noise severely impacts the results of using bandpass filtering. The energy of noise is notable in the time-frequency spectra (Figs 7c–f). The traditional and varimax norm based MR suppress a lot of noise, but damage signal significantly at the same time. In contrast, the RNMR method gets a successful result. Table 2 gives the quantitative analysis of the performances of the tested techniques. Figure 6. View largeDownload slide The second synthetic example. From left to right, traces 1–3: signal, Gaussian noise and limited-band random noise; trace 4: signal+Gaussian noise+limited-band random noise (input data); traces 5–11: seven scale components; traces 12–18: results using four bandpass filters (1–10–190–200 Hz, 25–35–165–175 Hz, 50–60–140–150 Hz, 75–85–115–125 Hz), MR, varimax norm based MR, and RNMR, respectively; traces 19–25: the corresponding removed noise (differences between trace 4 and traces 12–18); traces 26–32: the corresponding errors (differences between trace 1 and traces 12–18). Figure 6. View largeDownload slide The second synthetic example. From left to right, traces 1–3: signal, Gaussian noise and limited-band random noise; trace 4: signal+Gaussian noise+limited-band random noise (input data); traces 5–11: seven scale components; traces 12–18: results using four bandpass filters (1–10–190–200 Hz, 25–35–165–175 Hz, 50–60–140–150 Hz, 75–85–115–125 Hz), MR, varimax norm based MR, and RNMR, respectively; traces 19–25: the corresponding removed noise (differences between trace 4 and traces 12–18); traces 26–32: the corresponding errors (differences between trace 1 and traces 12–18). Figure 7. View largeDownload slide Comparison of time–frequency spectra of (a and b) clean and noisy data, (c–f) results after using 1–10–190–200 Hz, 25–35–165–175 Hz, 50–60–140–150 Hz, and 75–85–115–125 Hz bandpass filters, respectively, (g–i) results after using MR, varimax norm based MR, and RNMR, respectively. Figure 7. View largeDownload slide Comparison of time–frequency spectra of (a and b) clean and noisy data, (c–f) results after using 1–10–190–200 Hz, 25–35–165–175 Hz, 50–60–140–150 Hz, and 75–85–115–125 Hz bandpass filters, respectively, (g–i) results after using MR, varimax norm based MR, and RNMR, respectively. Table 2. Quantitative comparison of different approaches for the second synthetic example. BP represents bandpass filtering and Varimax MR represents varimax norm based MR.   Input  BP (1–10–190–200 Hz)  BP (25–35–165–175 Hz)  BP (50–60–140–150 Hz)  SNR  −13.0206  −7.6367  −7.4927  −6.5316    BP (75-85-115-125 Hz)  MR  Varimax MR  RNMR  SNR  −2.6882  −4.3178  −1.8567  4.6906    Input  BP (1–10–190–200 Hz)  BP (25–35–165–175 Hz)  BP (50–60–140–150 Hz)  SNR  −13.0206  −7.6367  −7.4927  −6.5316    BP (75-85-115-125 Hz)  MR  Varimax MR  RNMR  SNR  −2.6882  −4.3178  −1.8567  4.6906  View Large Real data examples For the development of shale gas, the hydraulic fracturing is an effective method to improve the permeability of rocks. However, the hydraulic fracturing is generally not strong and the accompanying seismicity is weak. Therefore we need an effective technique to detect the weak seismicity, then find out the location of the seismicity and further describe the fracture (e.g. the strike, length and thickness). The technique presented in this paper can be applied as a pre-processing technique to detect the weak signal stimulated in the process of hydraulic fracturing. To demonstrate the performance of the proposed approach in practice, we apply it to a real microseismic data set. The data set is recorded by 12 downhole 3-C receivers at temporal sampling of 0.5 ms and spatial sampling of 30 m. Fig. 8 gives the 3-D visualization of the spatial position of the monitoring system. The depth ranges from 1621 m to 1951 m for the receivers and from 2200 m to 3000 m for the hydraulic fracture growth. The horizontal distance between fracturing and monitoring wells is approximately 620 m. The position of perforation is X(North)=276.934 m, Y(East)=281.22 m, Z(Depth)=2023.12 m. The frequency band ranges from 0 to 350 Hz for the observed microseismic signals and from 0 to 900 Hz for the ambient noise. According to well logging data, the velocities of P wave and S wave range between [4000,6000] and [2500,3500], respectively. We approximate geological structure between 2200 m and 3000 m as a model of nine horizontal layers. Figure 8. View largeDownload slide 3-D visualization of the spatial position of the monitoring system. The red crosses, the black curve line and the blue circle represent the positions of the receivers, the fracturing well and the perforation, respectively. Figure 8. View largeDownload slide 3-D visualization of the spatial position of the monitoring system. The red crosses, the black curve line and the blue circle represent the positions of the receivers, the fracturing well and the perforation, respectively. There are totally eight injection stages in this project. The data used in this study is produced in the eighth stage. In this stage, the total monitoring time is approximately 3 hr. 205 microseismic events have been identified. The predicted length and width of the fracture grid are 920 m and 450 m, respectively. The recording time in the eighth stage is the longest in the whole project. As is often the case of microseismic monitoring, the signal induced from hydraulic fracturing is very weak compared with the ambient noise. Fig. 9 presents a typical 1500 ms 3-C record. The 1st–12th and 13th–24th traces are the two horizontal components data (H1 and H2), and the 25th–36th traces are the vertical component data (V). We observe that the initial data are extremely noisy, especially in the H1 and H2 records. The events are difficult to follow. Fig. 10 shows the amplitude spectrum of the initial record. The frequency band of the perforating signal in this stage ranges from 0 Hz to 500 Hz. Due to the impact of industrial electricity, there are low-frequency interferences in several traces. Figure 9. View largeDownload slide 3-C microseismic data. The 1st–12th traces correspond to the horizontal components (H1), the 13th–24th traces correspond to the horizontal components (H2), and the 25th–36th traces correspond to the vertical component (V). Figure 9. View largeDownload slide 3-C microseismic data. The 1st–12th traces correspond to the horizontal components (H1), the 13th–24th traces correspond to the horizontal components (H2), and the 25th–36th traces correspond to the vertical component (V). Figure 10. View largeDownload slide Amplitude spectra of initial (a) H1 component, (b) H2 component and (c) V component. Figure 10. View largeDownload slide Amplitude spectra of initial (a) H1 component, (b) H2 component and (c) V component. We first apply the bandpass filtering to the initial data. The results after using fixed one octave bandpass filtering, with trapezoidal bands 10–20–40–50 Hz, 30–40–80–90 Hz, 70–80–160–170 Hz, 150–160–320–330 Hz, and 310–320–640–650 Hz are demonstrated in Fig. 11, respectively. We can see that, the event in the V component is detectable in the 30–40–80–90 Hz (Fig. 11b) and 70–80–160–170 Hz (Fig. 11c) frequency bands, as the red arrows mark. The event in the H1 component is detectable in the 70–80–160–170 Hz (Fig. 11c) frequency band, as the red arrow marks. However, for the H2 component record, we can only find some faint coherent energy in the 70–80–160–170 Hz (Fig. 11c) frequency band, as the red arrow marks. But it is still difficult to identify the faint coherent energy as the microseismic signal. The bandpass filtering can improve the detectability of the weak signal in microseismic monitoring to some extent. Figure 11. View largeDownload slide Results after applying five bandpass filters. (a) 10–20–40–50 Hz component. (b) 30–40–80–90 Hz component. (c) 70–80–160–170 Hz component. (d) 150–160–320–330 Hz component. (e) 310–320–640–650 Hz component. Figure 11. View largeDownload slide Results after applying five bandpass filters. (a) 10–20–40–50 Hz component. (b) 30–40–80–90 Hz component. (c) 70–80–160–170 Hz component. (d) 150–160–320–330 Hz component. (e) 310–320–640–650 Hz component. We then decompose the initial data into five morphological scale components using various scales of semicircle SE, as shown in Fig. 12. As we observe, the coherent energy can be found more or less in the second to fourth scale components of both H1 and H2 component records, and in the second to fifth scale components of the V component. But we can hardly see any coherent energy in the 1st scale components of all H1, H2 and V component records and the fifth scale components of both H1 and H2 component records. Thus we choose the second to fourth scale components for the reconstructions of both H1 and H1 components, and second to fifth scale components for the reconstructions of the V component. We first compare the three different MR approaches. For the traditional MR, the weighting factor is chosen manually as (1, 1, 0.5), (1, 1, 0.5) and (1, 1, 1, 1) for H1, H2 and V components, respectively. Figs 13(a)–(c) demonstrate the three reconstructed results by RNMR, MR and varimax norm based MR. The RNMR approach obtains a satisfactory performance. We can easily follow the coherent energy in all H1, H2 and V components from Fig. 13(a). Both traditional and varimax norm based MRs suppress some ambient noise and make the events clearer than without processing, however, the proposed RNMR performs better. What if the traditional denoising techniques are applied to this data set? The next experiment is carried out by the well-developed MF (Justusson 1981) and SSA (Trickett 2008) methods. The MF approach is a well-known nonlinear approach which can smooth noise and enhance signal. The SSA is an up-to-date approach in seismics, which is implemented by using a rank-reduction operation, that is, the truncated singular value decomposition, to decompose the input data into an estimated noise subspace and an estimated signal subspace. The results after using MF and SSA are presented in Figs 13(d) and (e). We can see that, the median filtering performs well. The events are clearer than the initial data. But the SSA fails in the detection of weak signal, as we can see from Figs 13(d) and (e), in which events are still masked by the background noise. Fig. 14 presents the difference sections between the initial data (Fig. 9) and the processed results (Fig. 13). From the difference sections we can see that the RNMR, MR and MF approaches preserve the signal well. There is almost no coherent energy in Figs 14(a), (b) and (d). However, the Varimax norm based MR and SSA methods cause some signal leakages as indicated by the cyan arrows in Figs 14(c) and (e). Fig. 15 demonstrates a comparison of events picking in the initial data and results by our proposed RNMR, where the time picks are represented with red asterisks. It is clear that after using the proposed RNMR approach, the events become much clear and easy to pick. Fig. 16 shows a microseismic event locating result using the picked arrival times in V component after using RNMR, where the locating result is denoted by the cyan asterisk. The Geiger’s approach (Geiger 1912) is used to obtain the location. The calculation of travel time is based on the principle of ray tracing. The distribution of the major fracture generated in this injection stages (the eighth stage) has overall dimension of about 900 m along the strike and direction of 25 degree west by north. This located event, which is at the position of X = 316 m, Y = 291 and Z = 1813 m, reveals potential small-scale fractures generated with east-north strike direction. Figure 12. View largeDownload slide Multiscale morphological decomposition. (a) First scale component. (b) Second scale component. (c) Third scale component. (d) Fourth scale component. (e) Fifth scale component. Figure 12. View largeDownload slide Multiscale morphological decomposition. (a) First scale component. (b) Second scale component. (c) Third scale component. (d) Fourth scale component. (e) Fifth scale component. Figure 13. View largeDownload slide Results after using (a) RNMR, (b) MR, (c) Varimax norm based MR, (d) MF and (e) SSA. Figure 13. View largeDownload slide Results after using (a) RNMR, (b) MR, (c) Varimax norm based MR, (d) MF and (e) SSA. Figure 14. View largeDownload slide Difference sections between Fig. 9 and (a) Fig. 13(a), (b) Fig. 13(b), (c) Fig. 13(c), (d) Fig. 13(d) and (e) Fig. 13(e). Figure 14. View largeDownload slide Difference sections between Fig. 9 and (a) Fig. 13(a), (b) Fig. 13(b), (c) Fig. 13(c), (d) Fig. 13(d) and (e) Fig. 13(e). Figure 15. View largeDownload slide Arrival picking of (a) the initial data, and (b) the result after RNMR. Figure 15. View largeDownload slide Arrival picking of (a) the initial data, and (b) the result after RNMR. Figure 16. View largeDownload slide The results of locating. The cyan asterisk denotes the location of the microseismic event recorded in V components. Figure 16. View largeDownload slide The results of locating. The cyan asterisk denotes the location of the microseismic event recorded in V components. The proposed RNMR approach can be used not only in microseismic data, but also in reflection seismic data to detect weak signal. The complex subsurface structure, long travel time and offset always lead to a poor quality of data, where the reflected signal is with a low SNR and hardly to detect directly. Fig. 17 presents the deep part (3300–6700 ms) of a field shot gather. Due to the long travel time, the energy of the reflected signals is attenuated a lot. We can hardly follow the reflection events from the original profile (Fig. 17), except the two blurry events indicated by the two blue arrows. We apply the proposed and traditional MR as well as the two denoising tools, namely the SSA and MF approaches, trying to improve the quality of the original data and detect the weak signal. The processed results are shown in Fig. 18. Fig. 19 demonstrates the difference sections between the original data (Fig. 17) and processed results (Fig. 18). The proposed RNMR method gives us a very successful result. As we can see from Fig. 18(a), the two events are cleaner and smoother than those before reconstruction, as highlighted by the blue arrows. Moreover, we can also observe some weak events as highlighted by the red arrows, which are almost invisible in the original data. ALL the MR, MF and SSA approaches remove some background noise and improve the SNR of the original data. However, the weak events are still hardly to detect and follow, as shown in Figs 18(b)–(d). Figure 17. View largeDownload slide The deep part of a field shot gather. Figure 17. View largeDownload slide The deep part of a field shot gather. Figure 18. View largeDownload slide Results after using (a) RNMR, (b) MR, (c) MF and (d) SSA. Figure 18. View largeDownload slide Results after using (a) RNMR, (b) MR, (c) MF and (d) SSA. Figure 19. View largeDownload slide Difference sections between Fig. 17 and (a) Fig. 18(a), (b) Fig. 18(b), (c) Fig. 18(c) and (d) Fig. 18(d). Figure 19. View largeDownload slide Difference sections between Fig. 17 and (a) Fig. 18(a), (b) Fig. 18(b), (c) Fig. 18(c) and (d) Fig. 18(d). DISCUSSIONS In real microseismic and seismic data, a lot of factors, for example, complex subsurface structure, weak seismicity, long travel time and offset, may lead to weak signal. The proposed RNMR can be applied as a pre-processing technique to detect the weak signal and obtain its vital information such as the arrival time and amplitude. In this section, we will discuss the parameter selection, computational efficiency for RNMR in detail. Moreover, we will discuss how to extend RNMR from 1-D to 2-D or a higher dimensional version. All of these can help readers better implement the proposed RNMR method. Parameter selection There are two important parameters controlling the performance of RNMR, namely the morphological decomposition level and the smoothing radius of the shaping operator. The decomposition level controls the fineness of the multiscale morphological decomposition. Although the SE has three parameters: shape, height and width, all the differences in the three parameters will lead to the differences in morphological scales and ultimately result in the fineness of the multiscale decomposition. We adjust the three parameters (generally we fix the shape and only adjust height and width for one specific multiscale morphological decomposition) to change the morphological scale of the SE and then decompose the data into a number of multiscale components to obtain the morphological information. The larger the number of the decomposed components, the more careful the multiscale morphological analysis is. The performance of the RNMR method can be improved as the number of the decomposed components increases. However, the performance will converge to the optimum, when the decomposition number is large enough, for example, 7, 8 or 9 (because the fineness of the multiscale decomposition is enough to detect the weak signal). The smoothing radius controls the power of the shaping operator. The shaping operator poses a local-smoothness constraint to the weighting operator Wk (model) in RNMR. The greater the radius, the stronger the smoothing is. Fig. 20 shows a synthetic example, in which the RNMR approach with different input parameters is tested. Trace 1 is the clean data and trace 2 is the noisy data corrupted by Gaussian noise. Traces 3-10 are the results using different RNMRs with different input parameters, where, trace 3 corresponds to 3 scales and a radius of 9 samples; trace 4 corresponds to 5 scales and a radius of 9 samples; trace 5 corresponds to 7 scales and a radius of 9 samples; trace 6 corresponds to 9 scales and a radius of 9 samples; trace 7 corresponds to 11 scales and a radius of 9 samples; trace 8 corresponds to 7 scales and a radius of 3 samples; trace 9 corresponds to 7 scales and a radius of 15 samples; trace 10 corresponds to 7 scales and a radius of 21 samples. Traces 11-18 present the corresponding removed noise (differences between trace 2 and traces 3–10) and traces 19–26 present the corresponding error (differences between trace 1 and traces 3-10). Table 3 shows the corresponding quantitative analysis. Comparing traces 3-5, we can see that as the decomposition number increases, the performance of RNMR becomes better. However, there is no significant difference among traces 5–7, even though they have different decomposition numbers. Comparing traces 5 and 8–10, the RNMR approach is robust with respect to the choice of the smoothing radius, but may have a negative performance if we choose a radius that is too small (trace 8) or too large (trace 10). Figure 20. View largeDownload slide Comparison of the performance of RNMR with different parameters. From left to right, traces 1–2: clean and noisy data; traces 3–10: results using RNMR (3 scales, radius=9 samples), RNMR (5 scales, radius=9 samples), RNMR (7 scales, radius=9 samples), RNMR (9 scales, radius=9 samples), RNMR (11 scales, radius=9 samples), RNMR (7 scales, radius=3 samples), RNMR (7 scales, radius=15 samples), RNMR (7 scales, radius=21 samples), respectively; traces 11–18: the corresponding removed noise (differences between trace 2 and traces 3–10); traces 19–26: the corresponding error (differences between trace 1 and traces 3–10). Figure 20. View largeDownload slide Comparison of the performance of RNMR with different parameters. From left to right, traces 1–2: clean and noisy data; traces 3–10: results using RNMR (3 scales, radius=9 samples), RNMR (5 scales, radius=9 samples), RNMR (7 scales, radius=9 samples), RNMR (9 scales, radius=9 samples), RNMR (11 scales, radius=9 samples), RNMR (7 scales, radius=3 samples), RNMR (7 scales, radius=15 samples), RNMR (7 scales, radius=21 samples), respectively; traces 11–18: the corresponding removed noise (differences between trace 2 and traces 3–10); traces 19–26: the corresponding error (differences between trace 1 and traces 3–10). Table 3. Quantitative comparison of the performance of RNMR with different parameters. The units of radius are samples.   Input  RNMR (3 scales, radius=9)  RNMR (5 scales, radius=9)  SNR  −12.0489  4.7239  9.3527    RNMR (7 scales, radius=9)  RNMR (9 scales, radius=9)  RNMR (11 scales, radius=9)  SNR  14.2503  14.0604  13.9609    RNMR (7 scales, radius=3)  RNMR (7 scales, radius=15)  RNMR (7 scales, radius=21)  SNR  5.9451  12.6312  9.8808    Input  RNMR (3 scales, radius=9)  RNMR (5 scales, radius=9)  SNR  −12.0489  4.7239  9.3527    RNMR (7 scales, radius=9)  RNMR (9 scales, radius=9)  RNMR (11 scales, radius=9)  SNR  14.2503  14.0604  13.9609    RNMR (7 scales, radius=3)  RNMR (7 scales, radius=15)  RNMR (7 scales, radius=21)  SNR  5.9451  12.6312  9.8808  View Large Computational efficiency The computational issue of the presented RNMR will be discussed in the section. We test three parameters, namely the size of input trace, the smoothing radius and the number of decomposed components, that may influence the computation time of RNMR. The test results are shown in Fig. 21. Fig. 21(a) shows the computational costs of RNMR varying with the size of input trace, where the smoothing radius and the number of decomposed components are fixed to be 9 samples and 7, respectively. Fig. 21(b) shows the computational costs of RNMR varying with the smoothing radius, where the size of input trace and the number of decomposed components are fixed to be 1000 samples and 7, respectively. Fig. 21(c) shows the computational costs of RNMR varying with the number of decomposed components, where the size of input trace and the smoothing radius are fixed to be 1000 samples and 9 samples, respectively. It can be demonstrated that the computation time of RNMR increases exponentially as the size of input traces increases, and increases linearly as the number of decomposed components increases. The smoothing radius has little influence on the computational costs of RNMR. Figure 21. View largeDownload slide Computational costs analysis of RNMR. (a) Computing time costs vary with the size of input trace. (b) Computing time costs vary with the smoothing radius. (c) Computing time costs vary with the number of decomposed components. Figure 21. View largeDownload slide Computational costs analysis of RNMR. (a) Computing time costs vary with the size of input trace. (b) Computing time costs vary with the smoothing radius. (c) Computing time costs vary with the number of decomposed components. Extension of RNMR from 1-D to a higher dimensional version Essentially, the RNMR approach presented in this paper is a single-channel technique, which is applied on individual traces. However, a 2-D or higher dimensional version of RNMR is more powerful when two or more physical dimensions of the seismic data (to be processed) are available, because it can take more information of the input data into consideration. Here, we provide a potential way of extending RNMR from 1-D to a higher dimensional version. We use the extension of RNMR from 1-D to 2-D as an example. The other dimensional version of RNMR can be constructed in the similar way. The RNMR method can be extended from two aspects, a) from 1-D SE to 2-D SE and b) from 1-D shaping to 2-D shaping. The 2-D SE (or shaping) operator can be constructed by a tensor product of two 1-D SE (or shaping) operators with orthogonal directions as:   \begin{equation} {\bf b}_{xy}={\bf b}_{x}{\bf b}_{y}, \end{equation} (28)  \begin{equation} {\bf H}_{xy}={\bf H}_{x}{\bf H}_{y}, \end{equation} (29)where bxy and Hxy represent the 2-D SE and 2-D shaping operators, respectively. bx and by represent two 1-D SE operators with orthogonal directions. Hx and Hy represent two 1-D shaping operators with orthogonal directions. Generally speaking, the two orthogonal directions can be vertical (temporal) and horizontal (spatial) directions. Fig. 22 presents a simple synthetic model consisting of five traces. Fig. 22(a) is the clean data and Fig. 22(b) is the noisy data which is severely contaminated by Gaussian and limited-band random noise. We use this model as an example to show the different performance of RNMR working with different (SE and shaping) operators. Fig. 23 shows the comparison of the 3rd trace after using different MR based approaches. From left to right, traces 1–2 show clean and noisy data; traces 3–7 show results using traditional MR (working with 1-D SE), RNMR (working with 1-D SE and 1-D shaping), RNMR (working with 1-D SE and 2-D shaping), RNMR (working with 2-D SE and 1-D shaping), and RNMR (working with 2-D SE and 2-D shaping), respectively; traces 8–12 show the corresponding removed noise (differences between trace 2 and traces 3–7); traces 13–17 show the corresponding error (differences between trace 1 and traces 3–7). Table 4 provides the corresponding quantitative analysis. From this example, we can clearly see that different performance of the RNMR algorithm using or not using the spatial information. It can be demonstrated that the RNMR algorithm working with 2-D SE and 2-D shaping obtains the best result (trace 7). Figure 22. View largeDownload slide 2D synthetic example. Figure 22. View largeDownload slide 2D synthetic example. Figure 23. View largeDownload slide Comparison of the third trace of the synthetic 2-D data after using different approaches. From left to right, traces 1–2: clean and noisy data; traces 3–7: results using MR (1-D SE), RNMR (1-D SE, 1-D shaping), RNMR (1-D SE, 2-D shaping), RNMR (2-D SE, 1-D shaping), and RNMR (2-D SE, 2-D shaping), respectively; traces 8–12: the corresponding removed noise (differences between trace 2 and traces 3–7); traces 13–17: the corresponding error (differences between trace 1 and traces 3–7). Figure 23. View largeDownload slide Comparison of the third trace of the synthetic 2-D data after using different approaches. From left to right, traces 1–2: clean and noisy data; traces 3–7: results using MR (1-D SE), RNMR (1-D SE, 1-D shaping), RNMR (1-D SE, 2-D shaping), RNMR (2-D SE, 1-D shaping), and RNMR (2-D SE, 2-D shaping), respectively; traces 8–12: the corresponding removed noise (differences between trace 2 and traces 3–7); traces 13–17: the corresponding error (differences between trace 1 and traces 3–7). Table 4. Quantitative comparison of different approaches.   Input  MR (1-D SE)  RNMR (1-D SE, 1-D shaping)  SNR  −13.0206  −4.3178  4.6906    RNMR (1-D SE, 2-D shaping )  RNMR (2-D SE, 1-D shaping)  RNMR (2-D SE, 2-D shaping)  SNR  8.8904  3.3457  17.3369    Input  MR (1-D SE)  RNMR (1-D SE, 1-D shaping)  SNR  −13.0206  −4.3178  4.6906    RNMR (1-D SE, 2-D shaping )  RNMR (2-D SE, 1-D shaping)  RNMR (2-D SE, 2-D shaping)  SNR  8.8904  3.3457  17.3369  View Large CONCLUSIONS We have proposed a novel weak signal detection method based on mathematical morphological decomposition and reconstruction. We introduce a non-stationary weighting operator into the process of reconstruction, which can impel the reconstruction of weak signal. The RNMR is formulated as an inverse problem, which can be solved via the shaping regularization method. The proposed method is a decomposition-and-reconstruction technique. We first decompose the input data into a number of multiscale morphological components, and then reconstruct the data using the selected components incorporating with the non-stationary weighting operator. Removal of noise and detection of signal are simultaneously achieved when the reconstruction step is conducted. The detailed algorithm framework can help readers implement and further verify the presented theory. Synthetic and real data examples demonstrate its superior performance compared with the competing alternative approaches, that is, the bandpass filtering, median filtering, and SSA methods. Compared with the traditional MR method, an additional advantage of the proposed approach is that, it can adaptively choose the weighting coefficients for the reconstruction. Acknowledgements We would like to thank Yimin Yuan for the constructive discussion. We especially want to thank Dr. Ludovic Métivier and two anonymous reviewers for constructive suggestions that greatly improved the manuscript. YC would like to thank Sergey Fomel for inspiring discussions on the topic of shaping regularization. This work is supported by the National Basic Research Program of China (973 Program), grant NO: 2013CB228602. WH would like to thank the China Scholarship Council for the financial support. REFERENCES Bai M., Wu J., 2017. Efficient deblending using median filtering without correct normal moveout-with comparison on migrated images, J. Seism. Explor. , 26, 455– 479. Bolton H., Masters G., 2001. 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, J. geophys. Res. , 106( B7), 13 527– 13 540. https://doi.org/10.1029/2000JB900378 Google Scholar CrossRef Search ADS   Chen G., Fu L.Y., Chen K.F., Sun W., Wei W., Guan X., 2017. Calculation of the seismic imaging complexity of complex geological structures, J. Seism. Explor. , 26( 1), 81– 104. Chen H., Zhou H., Chen Y., Xie C., 2016a. A constant fractional-order Laplacian viscoacoustic wave equation, in 86th Annual International Meeting, SEG, Expanded Abstracts , pp. 3982– 3986. Chen W., Chen Y., Liu W., Cheng Z., 2016b. Nonstationary signal decomposition via improved complete ensemble empirical mode decomposition and its application in ground roll noise attenuation, in 86th Annual International Meeting, SEG, Expanded Abstracts , pp. 4139– 4144. Chen W., Chen Y., Xie J., Zu S., Zhang Y., 2016c. Multiples attenuation using trace randomization and empirical mode decomposition, in 86th Annual International Meeting, SEG, Expanded Abstracts , pp. 4498– 4502. Chen X., Wang R., Li H., Lu C., 2015. Application of multi-scaled morphology in microseismic weak signal detection, in 77th EAGE Conference and Exhibition 2015 , doi:10.3997/2214-4609.201413013. Chen Y., 2016. Dip-separated structural filtering using seislet thresholding and adaptive empirical mode decomposition based dip filter, Geophys. J. Int. , 206, 457– 469. https://doi.org/10.1093/gji/ggw165 Google Scholar CrossRef Search ADS   Chen Y., 2017. Fast dictionary learning for noise attenuation of multidimensional seismic data, Geophys. J. Int. , 209, 21– 31. Chen Y., 2018. Automatic microseismic event picking via unsupervised machine learning, Geophys. J. Int. , 212, 88– 102. https://doi.org/10.1093/gji/ggx420 Google Scholar CrossRef Search ADS   Chen Y., Fomel S., 2015a. Random noise attenuation using local signal-and-noise orthogonalization, Geophysics , 80( 6), WD1– WD9. https://doi.org/10.1190/geo2014-0227.1 Google Scholar CrossRef Search ADS   Chen Y., Fomel S., 2015b. EMD-seislet transform, in 85th Annual International Meeting, SEG, Expanded Abstracts , pp. 4775– 4778. Chen Y., Fomel S., Hu J., 2014. Iterative deblending of simultaneous-source seismic data using seislet-domain shaping regularization, Geophysics , 79, V179– V189. https://doi.org/10.1190/geo2013-0449.1 Google Scholar CrossRef Search ADS   Chen Y., Zhang D., Huang W., Chen W., 2016d. An open-source Matlab code package for improved rank-reduction 3D seismic data denoising and reconstruction, Comput. Geosci. , 95, 59– 66. https://doi.org/10.1016/j.cageo.2016.06.017 Google Scholar CrossRef Search ADS   Chen Y., Zhang D., Jin Z., Chen X., Zu S., Huang W., Gan S., 2016e. Simultaneous denoising and reconstruction of 5-D seismic data via damped rank-reduction method, Geophys. J. Int. , 206( 3), 1695– 1717. https://doi.org/10.1093/gji/ggw230 Google Scholar CrossRef Search ADS   Cheng Z., Chen W., Chen Y., Liu Y., Liu W., Li H., 2016. Improving the S-transform for high-resolution reservoir prediction and paleochannels delineation, in 86th Annual International Meeting, SEG, Expanded Abstracts , pp. 1742– 1747. Chouet B. et al.  , 2003. Source mechanisms of explosions at Stromboli volcano, Italy, determined from moment-tensor inversions of very-long-period data, in J. geophys. Res. , 108( B1), ESE 7-1–ESE 7-25. https://doi.org/10.1029/2003JB002535 Claerbout J.F., 1992. Earth Soundings Analysis: Processing versus Inversion , 6, Blackwell Scientific Publications. Ebrahimi S., Kahoo A.K., Chen Y., Porsani M.J., Chen W., 2016. A high-resolution weighted semblance for dealing with AVO phenomenon, in 78th Annual International Conference and Exhibition, EAGE, Extended Abstracts , doi:10.3997/2214-4609.201601327. Fehler M., House L., Phillips W.S., Potter R., 1998. A method to allow temporal variation of velocity in travel-time tomography using microearthquakes induced during hydraulic fracturing, Tectonophysics , 289( 1), 189– 201. https://doi.org/10.1016/S0040-1951(97)00315-6 Google Scholar CrossRef Search ADS   Fischer T., Hainzl S., Eisner L., Shapiro S., Le Calvez J., 2008. Microseismic signatures of hydraulic fracture growth in sediment formations: observations and modeling, J. geophys. Res. , 113( B2), doi:10.1029/2007JB005070. Flewelling S.A., Tymchak M.P., Warpinski N., 2013. Hydraulic fracture height limits and fault interactions in tight oil and gas formations, Geophys. Res. Lett. , 40( 14), 3602– 3606. https://doi.org/10.1002/grl.50707 Google Scholar CrossRef Search ADS   Fomel S., 2007a. Local seismic attributes, Geophysics , 72( 3), A29– A33. https://doi.org/10.1190/1.2437573 Google Scholar CrossRef Search ADS   Fomel S., 2007b. Shaping regularization in geophysical-estimation problems, Geophysics , 72( 2), R29– R36. https://doi.org/10.1190/1.2433716 Google Scholar CrossRef Search ADS   Fomel S., 2008. Adaptive multiple subtraction using regularized nonstationary regression, Geophysics , 74( 1), V25– V33. https://doi.org/10.1190/1.3043447 Google Scholar CrossRef Search ADS   Gan S., Wang S., Chen Y., Chen X., 2015a. Deblending of distance separated simultaneous-source data using seislet frames in the shot domain, in 85th Annual International Meeting , SEG Expanded Abstracts , pp. 65– 70. Gan S., Wang S., Chen Y., Chen X., 2015b. Seismic data reconstruction via fast projection onto convex sets in the seislet transform domain, 85th Annual International Meeting , SEG Expanded Abstracts , pp. 3814– 3819. Geiger L., 1912. Probability method for the determination of earthquake epicenters from the arrival time only, Bull. St. Louis Univ. , 8( 1), 56– 71. Gibbons S.J., Ringdal F., 2006a. The detection of low magnitude seismic events using array-based waveform correlation, Geophys. J. Int. , 165( 1), 149– 166. https://doi.org/10.1111/j.1365-246X.2006.02865.x Google Scholar CrossRef Search ADS   Gibbons S.J., Ringdal F., 2006b. The detection of low magnitude seismic events using array-based waveform correlation, Geophys. J. Int. , 165( 1), 149– 166. https://doi.org/10.1111/j.1365-246X.2006.02865.x Google Scholar CrossRef Search ADS   Gibbons S.J., Ringdal F., Kvaerna T., 2008. Detection and characterization of seismic phases using continuous spectral estimation on incoherent and partially coherent arrays, Geophys. J. Int. , 172( 1), 405– 421. https://doi.org/10.1111/j.1365-246X.2007.03650.x Google Scholar CrossRef Search ADS   Guo P., Mcmechan G.A., 2017. Evaluation of three first-order isotropic viscoelastic formulations based on the generalized standard linear solid, J. Seism. Explor. , 26( 3), 199– 226. Guo T., Chen H., Wang G., Wang H., Cheng M., Chen P., 2017. Seismic characteristics of complex halo-anhydrites and their effects on the underlying carbonatites: A case study of the right bank of the Amu Darya River, J. Seism. Explor. , 26( 4), 381– 397. Hajati T., Langenbruch C., Shapiro S., 2015. A statistical model for seismic hazard assessment of hydraulic-fracturing-induced seismicity, Geophys. Res. Lett. , 42( 24), 10 601– 10 606. https://doi.org/10.1002/2015GL066652 Google Scholar CrossRef Search ADS   Han J., Mirko V.D.B., 2015. Microseismic and seismic denoising via ensemble empirical mode decomposition and adaptive thresholding, Geophysics , 80( 6), KS69– KS80. https://doi.org/10.1190/geo2014-0423.1 Google Scholar CrossRef Search ADS   Hansen S.M., Schmandt B., 2015. Automated detection and location of microseismicity at Mount St. Helens with a large-N geophone array, Geophys. Res. Lett. , 42( 18), 7390– 7397. https://doi.org/10.1002/2015GL064848 Google Scholar CrossRef Search ADS   Hestenes M.R., Stiefel E., 1952. Methods of Conjugate Gradients for Solving Linear Systems , 49, NBS. House L., 1987. Locating microearthquakes induced by hydraulic fracturing in crystalline rock, Geophys. Res. Lett. , 14( 9), 919– 921. https://doi.org/10.1029/GL014i009p00919 Google Scholar CrossRef Search ADS   Huang W., Wang R., Chen Y., Yuan Y., Zhou Y., 2016a. Randomized-order multichannel singular spectrum analysis for simultaneously attenuating random and coherent noise, in 86th Annual International Meeting, SEG, Expanded Abstracts , pp. 4777– 4781. Huang W., Wang R., Zhou Y., Chen Y., Yang R., 2016b. Improved principal component analysis for 3D seismic data simultaneous reconstruction and denoising, in 86th Annual International Meeting, SEG, Expanded Abstracts , pp. 4777– 4781. Huang W., Wang R., Chen X., Chen Y., 2017a. Double least squares projections method for signal estimation, IEEE Transactions on Geoscience and Remote Sensing , 55( 7), 4111– 4129. https://doi.org/10.1109/TGRS.2017.2688420 Google Scholar CrossRef Search ADS   Huang W., Wang R., Chen X., Zhou Y., Chen Y., You J., 2017b. Low-frequency noise attenuation of seismic data using mathematical morphological filtering, in 87th Annual International Meeting, SEG, Expanded Abstracts , pp. 5011– 5016. Huang W., Wang R., Chen X., Zu S., Chen Y., 2017c. Damped sparse representation for seismic random noise attenuation, in 87th Annual International Meeting, SEG, Expanded Abstracts , pp. 5079– 5084. Huang W., Wang R., Zhang D., Zhou Y., Yang W., Chen Y., 2017d. Mathematical morphological filtering for linear noise attenuation of seismic data, Geophysics , 82( 6), 1– 78. https://doi.org/10.1190/geo2017-0321-TIOgeo.1 Google Scholar CrossRef Search ADS   Jia R.S., Zhao T.B., Sun H.M., Yan X.H., 2015. Micro-seismic signal denoising method based on empirical mode decomposition and independent component analysis, Chinese Journal of Geophysics , 58( 3), 1013– 1023. https://doi.org/10.1002/cjg2.20187 Justusson B., 1981. Median filtering: statistical properties, in Two-Dimensional Digital Signal Processing II , pp. 161– 196, ed. Huang T.S., Springer. Google Scholar CrossRef Search ADS   Kong S., Phinney R.A., Roy-Chowdhury K., 1985. A nonlinear signal detector for enhancement of noisy seismic record sections, Geophysics , 50( 4), 539– 550. https://doi.org/10.1190/1.1441931 Google Scholar CrossRef Search ADS   Koskinen L., Astola J.T., Neuvo Y.A., 1991. Soft morphological filters, in San Diego,’91, San Diego, CA , pp. 262– 270, eds Gader P.D., Dougherty E.R., International Society for Optics and Photonics. Li H., Wang R., Cao S., Chen Y., Huang W., 2016a. A method for low-frequency noise suppression based on mathematical morphology in microseismic monitoring, Geophysics , 81( 3), V159– V167. https://doi.org/10.1190/geo2015-0222.1 Google Scholar CrossRef Search ADS   Li H., Wang R., Cao S., Chen Y., Tian N., Chen X., 2016b. Weak signal detection using multiscale morphology in microseismic monitoring, J. appl. Geophys. , 133, 39– 49. https://doi.org/10.1016/j.jappgeo.2016.07.015 Google Scholar CrossRef Search ADS   Liu G., Fomel S., Chen X., 2011. Time-frequency analysis of seismic data using local attributes, Geophysics , 76( 6), P23– P34. https://doi.org/10.1190/geo2010-0185.1 Google Scholar CrossRef Search ADS   Liu W., Cao S., Chen Y., 2016. Application of variational mode decomposition in random noise attenuation and time-frequency analysis of seismic data, in 78th Annual International Conference and Exhibition, EAGE, Extended Abstracts , doi:10.3997/2214-4609.201601249. Lu B., Feng J., 2017. Coal working face imaging by seismic interferometry-using conveyer belt noise as source, J. Seism. Explor. , 26( 5), 411– 432. Maragos P., 1994. Morphological systems: slope transforms and max-min difference and differential equations, Signal Processing , 38( 1), 57– 77. https://doi.org/10.1016/0165-1684(94)90057-4 Google Scholar CrossRef Search ADS   Matheron G., 1975. Random Sets and Integral Geometry , John Wiley & Sons. Maxwell S., Rutledge J., Jones R., Fehler M., 2010. Petroleum reservoir characterization using downhole microseismic monitoring, Geophysics , 75( 5), 75A129– 75A137. https://doi.org/10.1190/1.3477966 Google Scholar CrossRef Search ADS   Michelet S., Toksöz M.N., 2007. Fracture mapping in the soultz-sous-forêts geothermal field using microearthquake locations, J. geophys. Res. , 112( B7), doi:10.1029/2006JB004442. https://doi.org/10.1029/2006JB004442 Mortazavi S.A., Javaherian A., Bidhendi M.N., Amindavar H.R., 2017. Application of adaptive pyramidal dual-tree directional filter bank in ground roll attenuation, J. Seism. Explor. , 26( 1), 49– 79. Mousavi S.M., Langston C.A., 2017. Automatic noise-removal/signal-removal based on general cross-validation thresholding in synchrosqueezed domain and its application on earthquake data, Geophysics , 82( 4), V211– V227. https://doi.org/10.1190/geo2016-0433.1 Google Scholar CrossRef Search ADS   Mousavi S.M., Horton S.P., Langston C.A., Samei B., 2016a. Seismic features and automatic discrimination of deep and shallow induced-microearthquakes using neural network and logistic regression, Geophys. J. Int. , 207, 29– 46. https://doi.org/10.1093/gji/ggw258 Google Scholar CrossRef Search ADS   Mousavi S.M., Langston C.A., Horton S.P., 2016b. Automatic microseismic denoising and onset detection using the synchrosqueezed continuous wavelet transform, Geophysics , 81( 4), V341– V355. https://doi.org/10.1190/geo2015-0598.1 Google Scholar CrossRef Search ADS   Neelamani R., Baumstein A.I., Gillard D.G., Hadidi M.T., Soroka W.L., 2008. Coherent and random noise attenuation using the curvelet transform, Leading Edge , 27( 2), 240– 248. https://doi.org/10.1190/1.2840373 Google Scholar CrossRef Search ADS   Pearson C., 1981. The relationship between microseismicity and high pore pressures during hydraulic stimulation experiments in low permeability granitic rocks, J. geophys. Res. , 86( B9), 7855– 7864. https://doi.org/10.1029/JB086iB09p07855 Google Scholar CrossRef Search ADS   Phillips W., Fairbanks T., Rutledge J., Anderson D., 1998. Induced microearthquake patterns and oil-producing fracture systems in the Austin Chalk, Tectonophysics , 289( 1), 153– 169. https://doi.org/10.1016/S0040-1951(97)00313-2 Google Scholar CrossRef Search ADS   Qu S., Verschuur D., Chen Y., 2016. Full waveform inversion using an automatic directional total variation constraint, in 78th Annual International Conference and Exhibition, EAGE, Extended Abstracts , doi:10.3997/2214-4609.201701340. Reshetnikov A., Buske S., Shapiro S., 2010. Seismic imaging using microseismic events: results from the San Andreas Fault System at SAFOD, J. geophys. Res. , 115( B12), doi:10.1029/2009JB007049. https://doi.org/10.1029/2009JB007049 Reyes-Montes J., Rietbrock A., Collins D., Young R., 2005. Relative location of excavation induced microseismicity at the underground research laboratory (AECL, Canada) using surveyed reference events, Geophys. Res. Lett. , 32( 5), doi:10.1029/2004GL021733. https://doi.org/10.1029/2004GL021733 Rutledge J.T., Phillips W.S., Schuessler B.K., 1998. Reservoir characterization using oil-production-induced microseismicity, Clinton County, Kentucky, Tectonophysics , 289( 1), 129– 152. https://doi.org/10.1016/S0040-1951(97)00312-0 Google Scholar CrossRef Search ADS   Rutledge J.T., Phillips W.S., Mayerhofer M., 2004. Faulting induced by forced fluid injection and fluid flow forced by faulting: an interpretation of hydraulic-fracture microseismicity, Carthage Cotton Valley Gas Field, Texas, Bull. seism. Soc. Am. , 94( 5), 1817– 1830. https://doi.org/10.1785/012003257 Google Scholar CrossRef Search ADS   Schimmel M., Paulssen H., 1997. Noise reduction and detection of weak, coherent signals through phase-weighted stacks, Geophys. J. Int. , 130( 2), 497– 505. https://doi.org/10.1111/j.1365-246X.1997.tb05664.x Google Scholar CrossRef Search ADS   Serra J., 1982. Image Analysis and Mathematical Morphology , 1, Academic Press. Shahin A., Myers M., Hathon L., 2017. Carbonates’ dual-physics modeling aimed at seismic reservoir characterization, J. Seism. Explor. , 26, 331– 349. Shapiro S., Dinske C., Rothert E., 2006. Hydraulic-fracturing controlled dynamics of microseismic clouds, Geophys. Res. Lett. , 33( 14), doi:10.1029/2006GL026365. Šílený J., Hill D.P., Eisner L., Cornet F.H., 2009. Non–double-couple mechanisms of microearthquakes induced by hydraulic fracturing, J. geophys. Res. , 114( B8), doi:10.1029/2008JB005987. Sinha D., Dougherty E.R., 1992. Fuzzy mathematical morphology, J. Vis. Commun. Image Represent. , 3( 3), 286– 302. https://doi.org/10.1016/1047-3203(92)90024-N Google Scholar CrossRef Search ADS   Sobania A., Evans J. P.O., 2005. Morphological corner detector using paired triangular structuring elements, Pattern Recognit. , 38( 7), 1087– 1098. https://doi.org/10.1016/j.patcog.2004.10.009 Google Scholar CrossRef Search ADS   Tian W., Wang Y., Zhou T., Xie C., Guan Z., Liu H., Xu K., 2017. Characteristics of high-frequency ultra-acoustic wave spectrum and pore size in low-permeability sandstone, J. Seism. Explor. , 26( 5), 399– 410. Tikhonov A., 1963. Solution of incorrectly formulated problems and the regularization method, Sov. Math. Dokl. , 5, 1035– 1038. Trickett S., 2008. F-xy cadzow noise suppression, in SEG Technical Program Expanded Abstracts 2008 , pp. 2586– 2590, Society of Exploration Geophysicists. Wang B., Wu R.-S., Chen X., Li J., 2015a. Simultaneous seismic data interpolation and denoising with a new adaptive method based on dreamlet transform, Geophys. J. Int. , 201( 2), 1180– 1192. Wang R., 2005. Noise-eliminated method by morphologic filtering in seismic data processing, Oil Geophysical Prospecting , 40( 3), 277– 282. Wang R., Wang Y., 2017. Seismic reflectivity inversion by curvelet deconvolution–a comparative study and further improvements , 26, pp. 331– 349. Wang R., Li Q., Zhang M., 2008. Application of multi-scaled morphology in denoising seismic data, Applied Geophysics , 5( 3), 197– 203. https://doi.org/10.1007/s11770-008-0033-3 Google Scholar CrossRef Search ADS   Wang Y., Zhou H., Li Q., Chen H., Gan S., Chen Y., 2015b. An unsplit convolutional perfectly matched layer for visco-acoustic wave equation with fractional time derivatives, in 85th Annual International Meeting , SEG Expanded Abstracts , pp. 3666– 3671, Society of Exploration Geophysicists. Wiggins R.A., 1978. Minimum entropy deconvolution, Geoexploration , 16( 1-2), 21– 35. https://doi.org/10.1016/0016-7142(78)90005-4 Google Scholar CrossRef Search ADS   Wu G., Fomel S., Chen Y., 2016. Data-driven time-frequency analysis of seismic data using regularized nonstationary autoregression, in 86th Annual International Meeting, SEG, Expanded Abstracts , pp. 1700– 1705. Xia Y., Ni S., Zeng X., 2013. Twin enigmatic microseismic sources in the gulf of guinea observed on intercontinental seismic stations, Geophys. J. Int. , 194, 362– 366. https://doi.org/10.1093/gji/ggt076 Google Scholar CrossRef Search ADS   Xie J., Di B., Wei J., Xie Q., Zu S., Chen Y., 2016. Stacking using truncated singular value decomposition and local similarity, in 78th Annual International Conference and Exhibition, EAGE, Extended Abstracts , doi:10.3997/2214-4609.201601325. Xie Q., Chen Y., Zhang G., Gan S., Wang E., 2015a. Seismic data analysis using synchrosqueezeing wavelet transform - a case study applied to Boonsville Field, in 77th Annual International Conference and Exhibition, EAGE, Extended Abstracts , doi:10.3997/2214-4609.201412752. Xie Q., Di B., Chen Y., Gan S., Jiao S., 2015b. Simultaneous structure enhancing filtering and fault detection using plane wave prediction, in 77th Annual International Conference and Exhibition, EAGE, Extended Abstracts , doi:10.3997/2214-4609.201412614. Yang W., Wang R., Chen Y., Wu J., 2014. Random noise attenuation using a new spectral decomposition method, 84th Annual International Meeting, SEG Expanded Abstracts , pp. 4366– 4370. Yuan Y., Wang R., Huang W., Chen X., Zhou Y., Jiang Y., 2016. Self-adaptive multi-scaled morphology for weak signal detection of thin interbedded reservoir, in 78th EAGE Conference and Exhibition 2016 . Zhang D., Chen Y., Gan S., 2016a. Iterative reconstruction of 3D seismic data via multiple constraints, in 78th Annual International Conference and Exhibition, EAGE, Extended Abstracts , doi:10.3997/2214–4609.201601241. Zhang D., Chen Y., Gan S., 2016b. Multidimensional seismic data reconstruction with multiple constraints, in 86th Annual International Meeting, SEG, Expanded Abstracts , pp. 4801– 4806. Zhang D., Chen Y., Huang W., Gan S., 2016c. Multi-step reconstruction of 3D seismic data via an improved MSSA algorithm, in CPS/SEG Beijing 2016 International Geophysical Conference & Exposition, SEG, Expanded Abstracts , pp. 745– 749. Zhang G., Wang Z., Sun Y., Chen Y., 2015. Characterization method of tight dolomitic reservoir of Fengcheng formation in the west slope of Mahu Sag, Junggar Basin, China, in 85th Annual International Meeting, SEG Extended Abstracts , pp. 2896– 2900. Zhang H., Zhang J., Li Z., Zhang J., Xiao J., 2017. Improving the imaging resolution of 3D PSTM in VTI media using optimal summation within the Fresnel zone, J. Seism. Explor. , 26( 4), 311– 330. Zhou Y., Wang R., Huang W., Yuan Y., Chen X., Wang L., Yang R., 2016. Attenuation of diffraction noise in marine surveys with mathematical morphology, in SEG Technical Program Expanded Abstracts 2016 , pp. 4644– 4648, Society of Exploration Geophysicists. Zu S., Zhou H., Chen Y., Liu Y., Qu S., 2015. A periodically variational dithering code for improving deblending, in 85th Annual International Meeting, SEG Extended Abstracts , pp. 38– 42. Zu S., Zhou H., Chen Y., Chen H., Cao M., Xie C., 2016a. A marine field trial for iterative deblending of simultaneous sources, 86th Annual International Meeting, SEG, Expanded Abstracts , pp. 113– 118. Zu S., Zhou H., Chen Y., Pan X., Gan S., Zhang D., Xie C., 2016b. Recovering the most from big gaps using least-squares inversion, 86th Annual International Meeting, SEG, Expanded Abstracts , pp. 4128– 4133. Zu S., Zhou H., Mao W., Zhang D., Li C., Pan X., Chen Y., 2017. Iterative deblending of simultaneous-source data using a coherency-pass shaping operator, Geophys. J. Int. , 211( 1), 541– 557. https://doi.org/10.1093/gji/ggx324 Google Scholar CrossRef Search ADS   © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society.
TI - Regularized non-stationary morphological reconstruction algorithm for weak signal detection in microseismic monitoring: methodology
JF - Geophysical Journal International
DO - 10.1093/gji/ggy054
DA - 2018-05-01
UR - https://www.deepdyve.com/lp/oxford-university-press/regularized-non-stationary-morphological-reconstruction-algorithm-for-vxZJAXxJwj
SP - 1189
EP - 1211
VL - 213
IS - 2
DP - DeepDyve
ER -