Unveiling the signals from extremely noisy microseismic data for high-resolution hydraulic fracturing monitoring

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.

techniques using different approaches such as: median filtering 23 , various kinds of mathematical transform based approaches [24][25][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.

Morphological Decomposition
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 are the morphological operations that process d with the SE τ b( ) as 41 : 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 χ d ) as: SCIENtIfIC REPORTS | 7: 11996 | DOI:10.1038/s41598-017-09711-2 We now use morphological opening and closing to represent data d.
, two indexed families of morphological opening and closing, respectively. Typically, the index k denotes the morphological scale. Whereupon, d can be represented as: Equations (5) and (6) can also be written as: 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 2nd-8th traces are the 7 scale components.

Traditional morphological reconstruction
For convenience, let: 1], are the morphological multi-scale components. The value of a sample t in c k is represented by ∈ c t f ( ) k . 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: where ∅ denotes empty set. The reconstruction of data by c k can be represented as: is the weighting coefficient that controls energy from different scale components. This decomposition allows for full reconstruction of the original data,

Orthogonalized morphological reconstruction
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: where ⋅ F 2 represents the squared Frobenius norm of a function. Equation (12) can be also represented as: . Thus, the orthogonalized morphological reconstruction holds as: 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 → c is on the line l. Operation Σ is a stretching transformation acting on vector → c . Equation (13) is actually to find the projection of vector → d on line l in the least-squares sense. Thus, an orthogonal decomposition of → d holds as: where ⋅ denotes Hadamard (or Schur) product. Hence, we name Σ orthogonalization operator. If we consider Σ → c as signal → s , and accordingly Σ → − → d c as background noise → n , equations (15) and (16) become the classical models used in [42][43][44][45] , Therefore, if we assume the weak signal is orthogonal to the background noise in microseismic monitoring, Σ → c is an estimation of the weak signal.

Solution of orthogonalization operator
The inversion problem in equation (13), however, is ill-posed. To stabilize the optimization, an extra regularization term is necessary to solve equation (13): where  represents the regularization operator. For convenience, we rewrite equation (19) as: 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: where ε is a scalar scaling parameter. The formal solution has the well-known form, where σ 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 : Combining equation (22) and (23), we have: By introducing scaling of A by λ 1/ in equation (24), we can rewrite it as: where λ is an introduced parameter controlling the physical dimensionality and enabling fast convergence when inversion is implemented iteratively 27 .

Implementation of the orthogonalized morphological reconstruction
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: 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: 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, Γ = HH T . Equation (25) can then be written as: The estimated weak signal s by the orthogonalized morphological reconstruction can be represented as:

Efficiency and effectiveness analysis of the orthogonalized morphological reconstruction
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.

Test of the orthogonalized morphological reconstruction
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 π/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 : 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 . (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 .  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.
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  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 (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-190 Hz, 70-80-120-130 Hz and 60-70-170-180 Hz, 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 . 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.3 ms and 0-200 Hz 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.

Application to a real data set
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.5 s record (8631-8632.5 s 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 14 5942 dB. 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  SCIENtIfIC REPORTS | 7: 11996 | DOI:10.1038/s41598-017-09711-2 ( 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 . 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: 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.
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.
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.

Conclusion
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     Table 1. Comparison of events detection in recorded data and denoised data after using different denoising approaches.
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.