Nothing Special   »   [go: up one dir, main page]

Next Article in Journal
Machine Learning and Hyperparameters Algorithms for Identifying Groundwater Aflaj Potential Mapping in Semi-Arid Ecosystems Using LiDAR, Sentinel-2, GIS Data, and Analysis
Next Article in Special Issue
Frequency-Domain Q-Compensated Reverse Time Migration Using a Stabilization Scheme
Previous Article in Journal
A Migratory Biomass Statistical Method Based on High-Resolution Fully Polarimetric Entomological Radar
Previous Article in Special Issue
A New Seismic Inversion Scheme Using Fluid Dispersion Attribute for Direct Gas Identification in Tight Sandstone Reservoirs
You seem to have javascript disabled. Please note that many of the page functionalities won't work as expected without javascript enabled.
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

An Automatic Velocity Analysis Method for Seismic Data-Containing Multiples

1
Key Laboratory of Geoghysical Exploration Equipment, Ministry of Education, Jilin University, Changchun 130026, China
2
College of Geo-Exploration Science and Technology, Jilin University, Changchun 130026, China
*
Author to whom correspondence should be addressed.
Remote Sens. 2022, 14(21), 5428; https://doi.org/10.3390/rs14215428
Submission received: 22 September 2022 / Revised: 26 October 2022 / Accepted: 27 October 2022 / Published: 28 October 2022
(This article belongs to the Special Issue Geophysical Data Processing in Remote Sensing Imagery)
Graphical abstract
">
Figure 1
<p>Comparison between seismic data without and with multiples and their velocity spectra, There is a corresponding relationship between event and peak marked with the same number in the figure, where the number 1 is primary, number 2 is surface-related multiple, and number is internal multiple: (<b>a</b>) seismic data containing only primaries. (<b>b</b>) Velocity spectra corresponding to data (<b>a</b>). (<b>c</b>) Seismic data containing primaries and multiples. (<b>d</b>) Velocity spectra corresponding to data (<b>c</b>).</p> ">
Figure 2
<p>A demonstration of the similarity between the original data velocity spectra and the predicted multiple velocity spectra: (<b>a</b>) the velocity spectra of original data. (<b>b</b>) The velocity spectra of predicted multiples. (<b>c</b>) Local similarity spectra.</p> ">
Figure 3
<p>Local similarity of the picked peaks.</p> ">
Figure 4
<p>Display of synthetic data: (<b>a</b>)Velocity model containing a velocity reversal, (<b>b</b>) CMP record of original data, (<b>c</b>) velocity spectra of original data.</p> ">
Figure 5
<p>Display of synthetic data test results: (<b>a</b>) muting-based strategy (MU), (<b>b</b>) Radon transform strategy (RT), (<b>c</b>) multiple attenuation strategy (MA), (<b>d</b>) proposed advanced method (AD).</p> ">
Figure 6
<p>Velocity comparison.</p> ">
Figure 7
<p>NMO comparison results. 1.4–2.4 s: (<b>a</b>) NMO results of the MP velocity, (<b>b</b>) NMO results of the AD velocity, (<b>c</b>) NMO results of the MA velocity. 3.2–4.2 s: (<b>d</b>) NMO results of the MP velocity, (<b>e</b>) NMO results of the AD velocity, (<b>f</b>) NMO results of the MA velocity.</p> ">
Figure 8
<p>Application on a single CMP gather from the field data: (<b>a</b>) CMP record of original data, (<b>b</b>) velocity spectra of original data, (<b>c</b>) velocity spectra and velocity curve of the control group, in which the red solid line is automatically obtained, and the yellow dotted line is the reference velocity, (<b>d</b>) velocity spectra and velocity curve obtained by the method proposed in this paper, in which the red solid line is automatically obtained, and the yellow dotted line is the reference velocity.</p> ">
Figure 9
<p>Comparison of velocity model: (<b>a</b>) reference model, (<b>b</b>) velocity model obtained by control group, (<b>c</b>) velocity model obtained by proposed method.</p> ">
Figure 10
<p>Comparison of stacked sections: (<b>a</b>) stack result of reference velocity, (<b>b</b>) stack result of the velocity obtained by control group, (<b>c</b>) stack result of the velocity obtained by proposed method.</p> ">
Review Reports Versions Notes

Abstract

:
Normal moveout (NMO)-based velocity analysis can provide macro velocity models for prestack data processing and seismic attribute inversion. Datasets with an increasing size require conventional velocity analysis to be transformed to a more automatic mode. The sensitivity to multiple reflections limits the wide application of automatic velocity analysis. Thus, we propose an automatic velocity analysis method for seismic data-containing multiples to overcome the limit of multiple interference. The core idea of the proposed algorithm is to utilize a multi-attribute analysis system to transform the multiple attenuation problem to a multiple identification problem. To solve the identification problem, we introduce the local similarity to attribute the predicted multiples and build a quantitative attribute called multiple similarity. Considering robustness and accuracy, we select two supplementary attributes based on velocity and amplitude difference, i.e., velocity variation with depth and amplitude level. Then we utilize the technique for order preference by similarity to ideal solution (TOPSIS) to balance weights for different attributes in automatic velocity analysis. An RGB system is adopted for multi-attributes fusion in velocity spectra for visualization and quality control. Using both synthetic and field examples to evaluate the effectiveness of the proposed method for data-containing multiples, the results demonstrate the excellent performance in the accuracy of the extracted velocity model.

Graphical Abstract">

Graphical Abstract

1. Introduction

Velocity is a critical parameter in both seismic data processing and imaging [1]. Methods have been proposed for velocity model-building such as full waveform inversion [2], joint-migration inversion [3], travel-time tomography [4], normal moveout based velocity analysis [5], and even machine learning [6]. Within these methods, normal moveout-based velocity analyses (hereinafter referred to as velocity analysis) are widely used in industry [7] due to their easy implementation, and no need for a prior velocity model. According to Wilson [8], there are two kinds of velocity analysis: one is to build velocity models for stacking the data for zero offset processing which can increase the signal-to-noise ratio of the stacked dataset [9,10]. The other is to obtain the moveout parameter model that can be used either to analyze the kinematic characteristics of the subsurface medium or as the initial model contributor for subsequent processes, such as travel time tomography, full waveform inversion [11], and time migration [12,13]. These applications make velocity analysis an indispensable step in seismic data processing [14]. Traditional velocity analysis requires processors to manually pick peaks in the velocity spectra. As shown in Figure 1a,b, the peak in the velocity spectra corresponds to the events in the seismic data. The coordinates of the peaks in the velocity spectra correspond to the kinematics (travel time and curvature) of the seismic wavefield, which means that the subsurface information can be obtained by picking manually the coordinate of the peaks in the velocity spectra. However, with the application of high-precision and even three-dimensional acquisition methods, in the face of massive data, manual picking in the velocity analysis can become a huge burden [15], which leads to that many automatic velocity analysis methods [16,17,18].
There are generally three ways of automatic velocity analysis. The first one is improving the quality of velocity spectra, which modifies the conventional semblance calculation [9] to increase the resolution of velocity spectra [19,20,21] and the adaptability of complex conditions [5,18]. The second one is optimizing the automatic picking algorithm and utilizing Monte Carlo [17], greedy algorithms [1], clustering analysis [22], simulated annealing [23] and other optimization algorithms to obtain high-precision automatic picking results. The third one is differential semblance velocity analysis (DSVA) [24,25], whose core conception is update migration velocity and flattening images automatically based on the residual moveout of neighboring traces in the image domain [26,27]. Because the accuracy of migration velocity depends on the approximation of the propagation theory to the seismic wavefield propagated in the underground medium, wave equation migration velocity analysis (WEMVA) has become popular [28,29].
Most automatic velocity analysis methods tend to chase higher peaks [30], which are sensitive to coherent noise, especially multiples [5,23,26]. Because the sea surface and the seabed are strong reflective interfaces, multiples cannot be avoided in marine seismic data. The kinematic characteristics of multiples are similar to those of primaries, so they appear as false peaks (Figure 1c,d) in the velocity spectra, which hinder the automatic velocity analysis methods from picking the true velocities. Compared with primaries, multiples reflect more than once between the sea surface and underground interface, and have more complex propagation paths. In Figure 1, two types of multiples are marked. The number two is marked as free surface multiples, and the number two is marked as an internal multiple. The existence of multiples will cause problems in automatic velocity analysis; for example, when using the automatic velocity analysis method to pick the false peaks generated by multiples, the obtained velocity model may be significantly different from the real formation velocity. Furthermore, when this unreliable model is used as the initial model of subsequent inversion such as FWI, it will not meet the accuracy requirements of the inversion [11], which will lead to inversion being trapped in a local non-informative minimum. So, the traditional seismic data processing requires multiple attenuations before velocity analysis. However, with the development of cognition and technology, in some industrial processes, multiples are regarded as effective signals to improve the imaging accuracy of complex structures, for example, imaging with multiples to reduce the relative amplitudes of the cross-talk noise [31] and enhance the amplitudes in the presence of strong scattering layers [32]. This means that velocity analysis of data containing multiples may become more common [33].
The existing automatic velocity analysis methods rely on two common strategies for dealing with multiples [12,34]. The first category is adopting a muting function in the velocity spectra based on the velocity difference, which assumes that multiples are low-velocity noise [35]. The muting-based strategy can be easily implemented [36], but not applicable to complex situations; for example, it can lead to small differences between primaries and multiples by the occurrence of shallow, high-velocity layers [37], which means the muting-based strategy is difficult to use to completely separate primaries and multiples [38]. The other method is to attenuate multiples before velocity analysis, through appropriate de-multiplying methods to obtain multiple-free input data. This approach can handle complex geological conditions and, thus has wider applicability [35,38]. The picking accuracy of the attenuation strategy depends on the success of multiple attenuation. Wave-equation-based multiple elimination (WEBME) is the most mainstream multiple attenuation method in the industry [39,40], which can be divided into two steps: multiple prediction and matching subtraction. Multiple prediction approximates the travel time of the different types of multiples, i.e., surface-related multiples [41], internal multiples [42], etc. Matching subtraction attenuates the prediction difference in amplitude and phase [43,44,45]. This means that obtaining multiple-free data is not a simple process, requiring a lot of manual intervention and computational costs [46]. Since currently published methods have some limitations in picking accuracy and automation, and multiples are unavoidable in marine exploration, studying an automatic velocity analysis method for seismic data containing multiples is meaningful.
In order to get rid of the requirement of automatic velocity analysis in primary-only input, we propose an automatic velocity analysis method for seismic data containing multiples. Firstly, we introduce the theory of traditional automatic velocity analysis and then describe how to reduce the sensitivity of automatic velocity analysis to multiples. The multi-attribute analysis system we built consists of multiple similarity, and velocity variation with depth and amplitude level to distinguish between primaries and multiples. Then, we adopt synthetic and field data examples to verify and evaluate the proposed method. Finally, we discuss the advantages and disadvantages of the proposed method, application recommendations and potential research directions.

2. Materials and Methods

2.1. Common Automatic Velocity Analysis

During automatic velocity analysis, the data-set ordered in shot gathers should be rearranged to common midpoint (CMP) gathers, and the velocity spectra corresponding to the data can be obtained through [9]
S ( v , t ) = j = t win t win i = 0 N 1 a [ k + j , x i ] 2 N j = t win t win i = 0 N 1 a [ k k + j , x i ] 2
where t and k is time index, k = t 2 + x i 2 v 2 , v is the scanning velocity, x i is the offset of trace i , a ( k + j , x i ) is the amplitude of the data, t win is the length of the time window, and N is the number of traces.
Then, the velocity curve L 1 [ v ( t ) ] can be automatically obtained by picking up the peaks in the velocity spectra s ( t , v ( t ) ) , for example, picking the peaks corresponding to the maximum of the variational integral [17]
L 1 [ v ( t ) ] = t min t max s ( t , v ( t ) ) d t ,
where t max , t min are the upper and lower boundaries of the time window, v ( t ) is the velocity corresponding to time sample t . As mentioned above, in the multiple attenuation process, it is necessary to make the predicted multiples match the actual situation in amplitude and phase to ensure that the multiple can be completely attenuated. Otherwise, we will need strict requirements in the selection of parameters and bounds to avoid the multiples, which is a non-adaptive and time-consuming task. Considering the calculation efficiency, we propose a method that can adaptively use the rough predicted multiples, which are obtained only by simple parameter adjustment and without matching subtraction, to perform automatic velocity analysis on seismic data containing multiples. To sum up, the key to this problem is how to eliminate the influence of multiple peaks on the automatic method by using the differential predicted multiples.

2.2. Multiple Independent Automatic Velocity Analysis

In order to eliminate the influence of multiples, we propose a multiple independent automatic velocity analysis method including the following steps:

2.2.1. Peak Picking

We first need to transform the data from shot gathers to velocity spectra. The original data containing multiples and the predicted multiples generated by modularization can be calculated by Equation (1) to obtain the corresponding velocity spectra. We denote the original data velocity spectra as D and the predicted multiple velocity spectra as M . As mentioned above, when the seismic data contain multiples, there will be many false peaks in the velocity spectra. Thus, it is unreasonable to pick higher peaks at each time step. In this paper, we choose the dichotomy method, pick all the peaks first, and then distinguish the peaks of primaries and multiples. Therefore, in this method we use a robust local maxima finder to find peaks of the velocity spectra [47]
φ = findpeaks ( D ) ,
where φ = { ( v , t ) | ( v 1 , t 1 ) , , ( v i , t i ) , , ( v z , t z ) } , i = 1 , 2 , , z and z is the number of peaks. The coordinate of each peak in the velocity spectra corresponds to its time sample indices and velocity, in which v , t expresses the velocity and time coordinates of peaks, respectively. This process can be easily implemented with the help of the FastPeakFind function in Matlab.

2.2.2. Attribution of Predicted Multiples

We found that, although the rough predicted multiple was obtained without fine parameter adjustment, the actual multiples in the original data are still similar to the predicted multiples, especially in the kinematic characteristics (travel time and curvature). The similarity is much greater than the similarity between the primary and the predicted multiples. In this case, the problem of how to distinguish primaries and multiples in the original data can be solved by measuring the similarity with the predicted multiples. To achieve this goal, this paper proposes a method based on local similarity to attribute the predicted multiples. We assume that the original data velocity spectra is composed of primaries, multiples and noise:
D = D p + D m + n 1 ,
where D p are primaries, D m are multiples, and n 1 is noise. Similarly, the predicted multiple velocity spectra can be written as:
M = M m + n 2 ,
where M m are the rough predicted multiples, n 2 is noise. Figure 2a,b show the velocity spectra of original data and the velocity spectra of the predicted multiple, where the red arrows indicate the primary locations and the purple arrows indicate the multiple location. We can see intuitively that the predicted multiples and the actual multiples match in number and locations, while the primaries do not. The reason for the above phenomenon is that the kinematic characteristics of the predicted multiples are similar to the actual multiples and orthogonal to the primaries. This difference is magnified in the velocity spectra.
In order to further analyze the similarity between the original data velocity spectra and predicted multiple velocity spectra, we introduce the local similarity for quantitative measurement [48]. Seismic data evaluated by local similarity are neither instantaneous nor global. Each datum has its corresponding local space; it can be summarized as a certain range, including an element and its adjacent elements [49]. The local similarity is given by:
s j = c 1 j c 2 j ,
c 1 j = δ [ λ 1 2 I + δ T A j T A j λ 1 2 I δ ] 1 δ T A j T b j ,
c 2 j = δ [ λ 2 2 I + δ T B j T B j λ 2 2 I δ ] 1 δ T B j T a j ,
where vector s j is local similarity between vectors a j and b j . a j and b j are the elements of the input data at trace index j respectively, j = 1 , 2 , , N , and N is the number of trace. c 1 j c 2 j denotes the Hadamard product between c 1 j and c 2 j . A j and B j are diagonal operators composed from the elements of a j and b j , that is, A j = diag ( a j ) , and B j = diag ( b j ) , δ is a triangle smoothing operator aiming to increase the smoothness in iterative optimization schemes [50,51], I is the identity matrix. λ 1 and λ 2 are two parameters controlling the physical dimensionality and enabling fast convergence when inversion is implemented iteratively. For a 2D case, the local similarity is calculated trace by trace, and finally arranged into a matrix S (Local similarity spectra), which has the same dimension as the input data, that is S = s 1 , s 2 , , s N .
Figure 2c shows the local similarity spectra between the predicted multiple and the original data. We can see that the local similarity highlights the similarity regions, in which high amplitude zones indicate multiples (purple arrows) and low amplitude zones indicate primaries (red arrows). In order to quantitatively show the similarity differences between multiples and primaries, we extracted the local similarity calculation results of all peak locations in Figure 2c, and the results are shown in Figure 3.
In Figure 3, the peaks are sorted according to their travel time. Among them, the first, second and fifth peaks are primaries, and the local similarities of these peaks are obviously different from the other peaks. Intuitively, the fourth peak and the seventh peak seem to have a low degree of similarity, but through the quantitative comparison, we can see that there is still an order of magnitude difference between the two peaks and the corresponding value of the primary peaks. In conclusion, even with the smallest multiple and largest primary, there is a difference of more than two orders of magnitude between them.
We calculated the maximum, minimum, average and median of their similarity (Table 1). They also prove that the primaries and multiples in the original data differ greatly in their similarity with the predicted multiple:
Note that, in principle, the predicted multiples do not include primaries, which means, compared with the primaries, multiples in the original data are much more similar to the predicted multiples. That is
L D m , M m L D p , M m ,
where L x , y denotes the local similarity between data x and y . To remove the multiple interference, we took this difference and introduced it into the automatic velocity analysis method. According to Equations (4) and (5)
L D , M = L D p + D m + n 1 , M m + n 2 .
We assume that the noise is distributed separately, then
L D , M = L D p + D m , M m = L D p , M m + L D m , M m .
According to Equation (9), Equation (11) can be further expressed as
L D , M = L D p , M m + L D m , M m L D m , M m ,
which means that the local similarity between the original data and the predicted multiples can be used as a reference to distinguish primaries without suppressing the multiples. Therefore, for any original data velocity spectra and predicted multiple velocity spectra, we can calculate their local similarity spectra to distinguish primaries and multiples; this equation is expressed as
S = L D , M
where S is the local similarity spectra between the predicted multiple and the original data. Considering the difference between the values of the primary and multiple in local similarity spectra is very large, we attribute them in the logarithmic form.
M S i = lg ( S ( t i , v i ) ) ,
M S i is the multiple similarity (MS) attribute of the i   th peak. lg ( S ( v i , t i ) ) = log 10 ( S ( v i , t i ) ) . Through Equation (12), we can quantitatively obtain the similarity between each peak and the predicted multiple.

2.2.3. Peaks Identification of Primary Reflection Based on the Multi-Attribute Analysis

The proposed multiple similarity provides a theoretical basis to select the peaks corresponding to the primary reflection from all the peaks. To improve the identification accuracy and robustness of the primary reflection peaks, we propose two supplementary attributes based on other multiple identification principles (velocity and amplitude), and introduce the multi-attribute analysis system [52] as a comprehensive evaluation tool to integrate different picking principles.
The multi-attribute analysis system can integrate and sort the criteria values of different strategies to achieve the selection of the suitable strategy under various qualitative/quantitative criteria [53]. This system mainly includes three steps: attribute definition and selection, determining attribute weights, and comprehensive evaluation [54]. To optimize the multi-attribute structure and reduce information redundancy as well as interference, it is necessary to ensure that all attributes have low correlation and relative independence in the multi-attribute analysis [55].
In this paper, we regard peaks of reflection wavefield (both primary and multiple) as different strategies of peaks, picking for the primary reflection and establishing a corresponding multi-attribute analysis system to decision-making in velocity analysis. The attributes used to distinguish between primary and multiple reflections in the velocity analysis mainly include velocity constraints [34], amplitude [56], bootstrapped differential semblance (BDS) [57] and detecting residual normal moveout [58]. Considering the calculation cost and the difficulty of adaptation, we selected the velocity and amplitude as supplementary attributes, and the three selected attributes are independent of each other, which ensures the stability of the multi-attribute analysis structure. The acquisition methods of the two supplementary attributes are as follows.
(1) Velocity variation with depth (VVD):
This attribute is mainly used to distinguish primaries and multiples through their velocity difference. The principle is that in general, the formation velocity increases with the increase in depth, while the propagation depth of multiples is less than that of primaries, which leads to the velocity of multiples being less than primaries. For any peak φ i , the velocity variation with depth, V V D i , can be expressed as
V V D i = v i V r e f t i ,
where v i , t i express are the velocity and time coordinates of peaks respectively, and V r e f t i represents the values corresponding to the velocity reference curve at the t i time sample indices. Furthermore, the velocity reference curve takes the coordinates of the velocity maximum peaks in different time windows as a reference and is obtained by interpolation along the time direction.
(2) Amplitude level (AL):
Amplitude difference is a classic attribute that mainstream automatic velocity analysis methods tend to use, as they pick the peaks with a higher amplitude [38]. Here, for peak φ i , its amplitude level attribute A L i is given by
A L i = lg D ( v i , t i ) lg max ( W i ) ,
where lg ( D ( v i , t i ) ) = log 10 ( D ( v i , t i ) ) , W i = { D ( v , t ) | t [ t i t w i n , t i + t w i n ] } , and t w i n represent the half width of the time window. Considering the thickness of the formation, the time window is set at 100 ms [59].
Note that the above three attributes have different units of measurement, value levels and physical meanings, which need to be regularized before determining attribute weights [60]. In multi-attribute analysis, different attributes can be divided into several categories, corresponding to different regularization equations. Among the three attributes applied by this method, VVD and AL are considered benefit attributes and MS is considered a cost attribute [61]. The regularization equations corresponding to the three attributes is as follows:
E L i = S M i S M min S M max S M min ,
E V i = V V D max V V D i V V D max V V D min ,
E A i = A L max A L i A L max A L min ,
where the values in the subscript indicate the maximum and minimum attribute values of all peaks, E L i indicates the value of the i   th peak under the condition of MS evaluation. Similarly, E V i and E A i indicate the value of the i   th peak under the condition of VVD and AL evaluation.
Through the above steps, we can regularize three different attributes. Then, this paper applies the technique for order preference by similarity to ideal solution (TOPSIS) method [62] to integrate the three attributes to automatically obtain the velocity curve. TOPSIS is convenient for application and calculation, and can output a quantitative evaluation [63]. The main steps are listed as follows [64]:
(1)
Determine the optimal peak φ o and the worst peak φ w ;
(2)
Obtain the relative proximity f ( i ) between each peak and the optimal peak;
(3)
Determine the dividing point of primary and multiple reflections through the inflection point of evaluation results [65].
In this method, there are z peaks to be evaluated and each peak has three evaluation criteria, whereby the smaller the value of each attribute after regularization, the closer to be primary the event is. Therefore, the optimal peaks φ o and the worst peaks φ w can be defined as
φ o = min E L 1 , E L 2 , E L z , min E V 1 , E V 2 , E V z , min E A 1 , E A 2 , E A z ,
φ w = max E L 1 , E L 2 , , E L z , max E V 1 , E V 2 , , E V z , max E A 1 , E A 2 , , E A z .
f ( i ) denotes the calculation of the relative proximity of the i   th peak to the optimal peak. It is calculated for each alternative and is defined as
f ( i ) = β i α i + β i i = 1 , 2 , , z
where α i represents the distance from the i   th peak to the optimal peak φ o , and β i represents the distance from the i   th peak to the worst peak φ w . We rank the peaks in descending order according to their corresponding f ( i ) . Through the equation, the primary peaks are those which are closest to the optimal peak and furthest from the worst peak [66]. f ( i ) is always between zero and one and one peak could be primary when it is closer to one [67]. α i and β i can be calculated by the following equations:
α i = ω L E L o E L i 2 + ω V E V o E V i 2 + ω A E A o E A i 2 ,
β i = ω L E L w E L i 2 + ω V E V w E V i 2 + ω A E A w E A i 2 ,
where E L o , E V o , E A o are the three attribute values of the optimal peak φ o , E L w , E V w , E A w are the three attribute values of the worst peak φ w , ω is the weight of the attribute. In this paper, we used the analytical hierarchy process (AHP) [68] method to determine the weight of each attribute. The specific choice is ω L = 0.6 , ω V = 0.2 , ω A = 0.2 , the detailed calculation process is highlighted in the Appendix A.
In addition, for conventional velocity analysis, quality control is difficult; once an error occurs, the processor needs to check all the peaks, which is a large burden. Therefore, this method visualizes the automatic velocity analysis process by using an RGB system [69] during processing. The three attributes of the peak are captured into a two-dimensional chart by color. Different colors in the RGB spectra represent different attributes. In this way, the difference between multiples and primaries can be directly displayed by color, so that the processor can focus on key peaks, which reduces the cost and difficulty of manual correction. The process of the proposed method shown in Algorithm 1.
Algorithm 1: Multiple Independent Automatic Velocity Analysis Algorithm
(1) Input: Original data containing multiples d , Predicted multiple generated by modularization m
(2) Calculate the velocity spectra for: D , M
(3) Peak picking in velocity spectra D : φ findpeaks ( D )
φ = { ( v , t ) | ( v 1 , t 1 ) , , ( v i , t i ) , , ( v z , t z ) }
(4) Attribute predicted multiples: S = L D , M
(5) i 1
(6) Main cycle:
(9)   Acquire MS attribute: M S i = lg ( S ( t i , v i ) )
(10)    Acquire VVD attribute: V V D i = v i V r e f t i
(11)    Acquire AL attribute: A L i = lg s d ( v i , t i ) lg max ( W i )
(12)     i i + 1
(13) Until i = z
(14) Regularization: E L i = S M i S M min S M max S M min , E V i = V V D max V V D i V V D max V V D min , E A i = A L max A L i A L max A L min
(15) Determine primary based on Multi-attribute Analysis theory:
Calculate relative proximity: f ( i ) = β i α i + β i , i = 1 , 2 , , z
Determine the dividing point between primaries and multiples (Inflection point): f ( i ) = 0
(16) RGB space mapping
(17) Output: Velocity curve, RGB velocity spectra

3. Results

A synthetic example was first used to test the proposed method. In order to increase the variety and complexity of multiples, we used a layered model containing a velocity reversal layer as shown in Figure 4a. The CMP gather is shown in Figure 4b, and Figure 4c is the corresponding velocity spectra. In the CMP gather, we can see that primaries and multiples share the same kinematic characteristics and strong amplitudes, making it difficult to distinguish primaries, especially for 2.5–4.5 s. In the velocity spectra shown in Figure 4c, due to the influence of the low-velocity layer (1.5–2.25 km), the velocities of multiples and primaries are close to each other. Thus, false peaks and primary peaks have similar velocities in the range of 2.5–4.5 s. These false peaks will largely interfere with automatic methods of picking velocities.
We applied four different automatic velocity analysis strategies in the test of synthetic data; the results are shown in Figure 5. Figure 5a is the result using muting-based strategies (we refer to it as MU), Figure 5b shows the result by radon transform [1] (we refer to it as RT), Figure 5c displays the result based on multiple attenuation strategies (we refer to it as MA), and Figure 5d is the result through the advanced multiple independent method (we refer to it as AD). By comparing the Figure 5 with Figure 4b, we can have a qualitative understanding of the multiple processing capabilities of various strategies. In Figure 5a, The white curve represents the muting range. It can be seen that in the shallow time zone (before 2.5 s), when multiples and primaries are scattered, the muting strategy can filter out the multiples, but in the range of 2.5–4.5 s, since multiple and primary are difficult to distinguish, most multiples are retained to avoid damaging the effective signal. For RT strategy in Figure 5b, the sparse hyperbolic Radon transform provides a sparse domain for velocity analysis; some multiples are attenuated in the transformation process, but there are still considerable residues (especially after 2.5 s). Comparing Figure 5c with Figure 4b, after multiple attenuation, most of the false peaks generated by multiples are attenuated, but there are still small amounts of multiple residuals. Figure 5d is the velocity spectra corresponding to proposed method, in which the peaks of primaries and multiples are displayed in color. In Figure 5, the red curves on each figure are the velocities obtained by the corresponding method. The greater the interference of multiples, the mor fluctuation the velocity curve will have. In order to evaluate the effectiveness of the proposed methods, we have made a quantitative analysis of the four velocity curves.
We manually picked the peaks of primaries referring to the true velocity model, and using the manual picking (refer to it as MP) result as the reference velocity. To test the effectiveness of the proposed method, we compared the velocities obtained by various strategies with the reference curve, and used σ to quantitatively measure the difference between the reference and estimated velocities, The results are shown in Figure 6 and Table 2; the equation of σ is as follows:
σ = v r v e 2 v r 2 ,
where v r is reference velocity and v e is estimated velocity. In Figure 6, the yellow, black, brown, blue and red lines represent MP, MU, RT, MA and AD velocities, respectively. Corresponding to the quantitative analysis results in Table 2, we can see that the estimated velocity by AD is the closest one to reference and that the red and yellow curves are almost overlapping each other. Correspondingly, Table 2 shows that between the AD velocity and the reference velocity, the calculated metric is σ = 1.3 % which is the smallest among the four strategies. The estimated velocity using the MA strategy is suboptimal, for which σ = 5.0 % . The main reason for the impact on the estimation accuracy can be clearly found in Figure 6 and Figure 5c, That is, the residual multiples at 2.5 s makes the automatic method pick a false peak. Additionally, more false picks are picked by the RT method and MU method, which further reduces the accuracy of their estimated velocities, for which the σ are 5.4% and 6.0% respectively. This comparison result shows the dependence of traditional automatic velocity analysis algorithm on multiple attenuations. However, the method proposed in this paper can measure the difference between the actual situation and the predicted multiples through local similarity, effectively improving the tolerance.
In order to test the accuracy of the velocity obtained by proposed method, we respectively used the reference velocity (MP), the velocity obtained by the proposed AD method, and the velocity obtained by the MA strategy, which is suboptimal to perform NMO correction on the original data. For a closeup look, we zoom the data in time regions 1.4–2.4 s and 3.2–4.2 s. The results are shown in Figure 7.
Within 3.2–4.2 s, the primary events were basically flattened (marked by yellow dotted line). However, in the 1.4–2.4 s, due to the influence of the residue false peak, the velocities obtained by MA strategy were lower than the real ones, which led to the overcorrection of the multiple events (marked by red dotted line). So, the results of NMO correction proved the effectiveness of the proposed method for data containing multiples.
In order to further verify the effectiveness of the proposed method, we applied it to field marine data for testing. The proposed method first was applied to a single CMP gather for quantitative analysis, and then we used the proposed method to process all CMP gathers to obtain the 2D velocity model of field data. For the data, the CMP interval is 25 m, the sampling interval is 4 ms, and maximum time is 2.8 s. Further, in the test on field data, because the real velocity model is unknown, we take the velocity curve picked manually as reference. As the multiple attenuation strategy can achieve better results than the other two strategies, we take the multiple attenuation strategy as the comparison object of the proposed method in the test of field data.
By comparing Figure 8b,c, we can see that after the multiple attenuation processing, the effective signals in the field data have been highlighted, which improved the data quality, but multiple residues can still be observed. The false peaks generated by these multiples will seriously interfere with the automatic velocity analysis. Comparing the velocity curves in Figure 8c,d, we can see that the velocity obtained by the proposed method is closer to the reference velocity. In order to further test the stability of the method, we applied it to the entire data to obtain a velocity model. The velocity model is shown in Figure 9, (a) is the reference model picked manually, (b) is the velocity model obtained by the control group, and (c) is the velocity model obtained by the proposed method. There is an obvious gap in the 1.2–1.6 s area (marked by red dotted line). The reason is that the model obtained by the control group is affected by residual multiples in this area. Meanwhile, in Figure 9b, there is a mixture of high-velocity and low-velocity formation (marked by red arrow), as a consequence of the incorrect picking of low-velocity multiple peaks by the control method.
We respectively applied the three velocity models to the original data to calculate their stacked sections. Figure 10a shows the stack section corrected with reference velocity, (b) is stack section of the control groups and (c) is the stack section with velocity obtained by the proposed method. In comparison, one event exists in Figure 9a,c near 1.6 s (marked by yellow dotted line), but in the stacked section of the control group, this event cannot be observed. The reason is that near 1.6 s, the control method is affected by the residual multiple and obtained the wrong velocity, so that the primary event was covered. In the comparison of three stacked sections, the stacking sections with velocities obtained by the proposed method are basically consistent with the reference-stacked section, which proves the effectiveness of the method for data-containing multiples.

4. Discussion

NMO-based velocity analysis is an indispensable step for whole seismic exploration, and can provide the initial velocity model for other velocity analysis methods based on the inversion problem [11,13]. The traditional interactive picking method may not meet the processing requirements of massive seismic data. However, the effectiveness of the published automatic methods is degraded by multiple interference. This gap for automatic velocity analysis makes the subject of the proposed method meaningful. Compared with conventional automatic velocity analysis, the proposed method significantly eliminates the influence of multiples and obtain an accurate velocity curve, which has been demonstrated in the synthetic and field examples.

4.1. Sensitivity to Multiple

For multiple-free datasets, conventional automatic velocity analysis can obtain a relatively accurate picked velocity in a fully automatic way. The coherent noise in CMPs greatly affects the peaks picking, leading to a lower picked velocity. To overcome this limit, several methods (e.g., multiple attenuation before velocity analysis [35], muting low-velocity zone in velocity spectra [38], interactive guidance [70], dip filter in CMPs [26]) are proposed, which make the current velocity analysis semi-automatic. In this paper, as to achieve velocity analysis for data-containing multiples in a fully automatic way, we attributed the predicted multiples and combined them with two additional auxiliary attributes through multi-attribute analysis as the novel principle of picking peaks.

4.2. Processing Detail

For velocity picking, the accuracy of picked velocity depends on the resolution of the velocity spectra. For improving the resolution, many researchers propose several modified semblances by measuring the difference between the neighbor traces in some domains and building the corresponding time-weighting function, such as AB semblance [18], offset-dependent weighting semblance [5], similarity-weighted semblance [14] and so on. For the proposed method, we adopt the traditional semblance proposed by Taner [9]. The reason for choosing the traditional one is automation and computing efficiency. The resolution of these modified semblance depends on the applicability of the weighting function, which means parameter adjustment and additional calculation costs. However, we still recommend utilizing an advanced semblance calculation in complex situations to improve the stability of automatic velocity analysis, such as data with the amplitude variation with offset (AVO) anomaly.
Among the attributes proposed in this paper, multiple similarity is the most critical. Its basic assumption is that the existing methods can accurately predict the travel time of the multiples. In order to meet the assumption as much as possible, we recommend predicting different types of multiples (e.g., surface-related multiples, internal multiples, water-bottom multiples) by traditional multiple prediction methods in the industry. In this paper, we highlight using the surface-related multiple elimination (SRME) approach [41] to predict surface-related multiples and common-focus-point (CFP) approach [71] to predict internal multiples, etc. The above strategies are obtained by balancing prediction and computational efficiency. In theory, the completeness of predicting multiples determines the accuracy of identifying multiples. We adopted multiple prediction without adjusting parameters to efficiently obtain different types of multiples. Then, we utilized local similarity to ensure the accuracy of the pro-posed automatic velocity analysis under the condition of large prediction error. It is true that when the formation conditions are complex, some special multiples such as diffraction multiples cannot be predicted. At this time, the multiple similarity attributes may not be applicable, but the other two auxiliary attributes can be combined to achieve primary extraction.

4.3. Usage Recommendation

For the proposed method, we recommend two types of uses:
(1)
Batch processing when the data size is too large. For instance, when processing three-dimensional data, automatic velocity analysis can reduce a lot of burden in the face of massive CMPs. The proposed method is fully automatic and multiple-independent, which is more suitable for adaptive batch processing.
(2)
Instead of the traditional semblance spectra, the RGB-based color spectra provide a more intuitive distinction between primaries and multiples. The color of each peak implies rich geophysical information, which can be used as a reference for manual picking. The relationship between seismic wavefield and color is shown in Table 3:
The RGB display allows a more convenient interaction mode in peaks picking. The processor can just focus on the suspicious peaks and obtain a more accurate velocity model with the assistance of color.

4.4. Potential Research Directions

Recently, deep learning has rapidly developed and even surpassed the human level in dealing with some practical problems, such as real-time processing. Automatic seismic inversion and processing with deep learning has gradually become a trend [6,72]. However, the training process requires large numbers of labeled samples [73]. For NMO-based velocity analysis, a standardized processing flow is needed, such as pre-processing, multiple attenuation and manual picking, which are time-consuming and laborious. A large number of high-quality training sample sets can be obtained through the proposed method combined with manual inspection.
Multiples are considered a contribution to imaging, which increases demand for velocity model building of dataset-containing multiples. The proposed method can match this demand and provide a micro-velocity model for follow-up imaging application. In this way, migration velocity analysis based on the wave equation can be extended to data-containing multiples.

5. Conclusions

We propose a multiple independent automatic velocity analysis method to solve the problem that automatic picking methods do not work with multiples. The main idea here is the utilization of the concept of multiple similarity as a novel peak-picking principle to break the primary dominant assumption. We combine multiple similarity with classical attributes such as velocity and amplitude as a compound-picking principle and adopt multi-attribute analysis theory to pick adaptively the peaks of the primaries. The proposed method has been applied to synthetic and field data tests, and has achieved better results than the control group, which verifies its effectiveness to data-containing multiples.
Another advantage of the proposed method is the visualization. Different from the traditional semblance spectra, we connect the three proposed attributes to the RGB system and display these peaks by color, leading to one peak that can demonstrate more comprehensive and intuitive geophysical information for processers.
The proposed method has the advantages of automation, efficiency and is easy to implement. Considering the future development trend of seismic exploration, the proposed multiple independent automatic velocity analysis method has practical value and sufficient development potential.

Author Contributions

All authors made significant contributions to this paper. J.Z.: algorithm writing, data analysis and original manuscript writing. B.H.: Investigation, development of ideas and review the manuscript. D.W.: Supervision, conceptualization, and field data acquisition. X.G.: Modeling, data test and manuscript editing. All authors have read and agreed to the published version of the manuscript.

Funding

This research was supported in part by the National Natural Science Foundation of China (No. 42074151, No. 41374108), and Major Projects of the National Science and Technology of China (Grant No. 2016ZX05026-002-003). This project is supported by special fund of Key Laboratory of Geoghysical Exploration Equipment, Ministry of Education (Jilin University).

Data Availability Statement

Not applicable.

Acknowledgments

Thanks for the open source program of Matlab.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

Analytic hierarchy process (AHP) [74] is used as the common weight determination method in multi-attribute Analysis to determine the weight of each evaluation attribute and the determination process is as shown in the equation
O = 1 o 1 o 2 1 o 1 1 o 3 1 o 2 1 o 3 1 ,
where O represents the decision matrix, o 1 is the scale of the first attribute relative to the second attribute, o 2 is the scale of the first attribute relative to the third attribute, and o 3 is the scale of the second attribute relative to the third attribute.
The specific value of the decision matrix can be obtained by looking up the scale table in the AHP method [74]. AHP method does not need to quantify the importance of each attribute but it only needs to qualitatively determine the importance of each attribute to obtain the scale of each attribute. The VVD and AL attributes are more likely to be affected by complex terrain because they are based on the difference of velocity and amplitude. Because the multiple prediction method is widely used, therefore, it is qualitatively believed that the MS attribute is slightly more important than the others. Moreover, the decision matrix can be converted according to the scale table of AHP method into the form shown in (A2).
O = 1 3 3 1 3 1 1 1 3 1 1 .
After the establishment of the decision matrix, the weight of each attribute in the AHP method can be obtain by normalizing the principal vector (eigen vector) of the matrix [73]. Through this way, it can be found that ω L = 0.6 , ω V = 0.2 , ω A = 0.2 .

References

  1. Chen, Y. Automatic velocity analysis using high-resolution hyperbolic Radon transform. Geophysics 2018, 83, A53–A57. [Google Scholar] [CrossRef]
  2. Virieux, J.; Asnaashari, A.; Brossier, R.; Métivier, L.; Ribodetti, A.; Zhou, W. An introduction to full waveform inversion. In Encyclopedia of Exploration Geophysics; Society of Exploration Geophysicists: Houston, TX, USA, 2017; pp. R1-1–R1-40. [Google Scholar] [CrossRef]
  3. Verschuur, D.; Staal, X.; Berkhout, A. Joint migration inversion: Simultaneous determination of velocity fields and depth images using all orders of scattering. Lead. Edge 2016, 35, 1037–1046. [Google Scholar] [CrossRef]
  4. Tanis, M.C.; Shah, H.; Watson, P.A.; Harrison, M.; Yang, S.; Lu, L.; Carvill, C. Diving-wave refraction tomography and reflection tomography for velocity model building. In Proceedings of the 2006 SEG Annual Meeting, New Orleans, LA, USA, 1–6 October 2006. [Google Scholar] [CrossRef] [Green Version]
  5. Luo, S.; Hale, D. Velocity analysis using weighted semblance. Geophysics 2012, 77, U15–U22. [Google Scholar] [CrossRef]
  6. Kazei, V.; Ovcharenko, O.; Plotnitskii, P.; Peter, D.; Zhang, X.; Alkhalifah, T. Mapping full seismic waveforms to vertical velocity profiles by deep learning. Geophysics 2021, 86, R711–R721. [Google Scholar] [CrossRef]
  7. Khoshnavaz, M.J. High-resolution seismic velocity analysis by sign-based weighted semblance. Geophysics 2021, 86, U135–U143. [Google Scholar] [CrossRef]
  8. Wilson, H.; Gross, L. Reflection-constrained 2D and 3D non-hyperbolic moveout analysis using particle swarm optimization. Geophys. Prospect. 2019, 67, 550–571. [Google Scholar] [CrossRef]
  9. Taner, M.T.; Koehler, F. Velocity spectra—Digital computer derivation applications of velocity functions. Geophysics 1969, 34, 859–881. [Google Scholar] [CrossRef]
  10. Alkhalifah, T. Seismic data processing in vertically inhomogeneous TI media. Geophysics 1997, 62, 662–675. [Google Scholar] [CrossRef]
  11. Virieux, J.; Operto, S. An overview of full-waveform inversion in exploration geophysics. Geophysics 2009, 74, WCC1–WCC26. [Google Scholar] [CrossRef]
  12. Yilmaz, Ö. Seismic Data Analysis: Processing, Inversion, and Interpretation of Seismic Data; Society of Exploration Geophysicists: Houston, TX, USA, 2001. [Google Scholar] [CrossRef]
  13. Zhou, H.-W.; Hu, H.; Zou, Z.; Wo, Y.; Youn, O. Reverse time migration: A prospect of seismic imaging methodology. Earth-Sci. Rev. 2018, 179, 207–227. [Google Scholar] [CrossRef]
  14. Chen, Y.; Liu, T.; Chen, X. Velocity analysis using similarity-weighted semblance. Geophysics 2015, 80, A75–A82. [Google Scholar] [CrossRef]
  15. Li, Y.; Zhao, D.; Yu, S. Status and trend on land seismic acquisition technique of SINOPEC. Geophys. Prospect. Pet. 2013, 52, 363–371. (In Chinese) [Google Scholar]
  16. Toldi, J.L. Velocity analysis without picking. Geophysics 1989, 54, 191–199. [Google Scholar] [CrossRef]
  17. Lumley, D.E. Monte Carlo automatic velocity picks. SEP Rep. 1997, 75, 1–25. [Google Scholar]
  18. Fomel, S. Velocity analysis using AB semblance. Geophys. Prospect. 2009, 57, 311–321. [Google Scholar] [CrossRef] [Green Version]
  19. Gan, S.; Wang, S.; Chen, Y.; Qu, S.; Zu, S. Velocity analysis of simultaneous-source data using high-resolution semblance-coping with the strong noise. Geophys. J. Int. 2016, 204, 768–779. [Google Scholar] [CrossRef] [Green Version]
  20. Gong, X.; Wang, S.; Zhang, T. Velocity analysis using high-resolution semblance based on sparse hyperbolic Radon transform. J. Appl. Geophys. 2016, 134, 146–152. [Google Scholar] [CrossRef]
  21. Ebrahimi, S.; Kahoo, A.R.; Chen, Y.; Porsani, M. A high-resolution weighted AB semblance for dealing with amplitude-variation-with-offset phenomenon. Geophysics 2017, 82, V85–V93. [Google Scholar] [CrossRef]
  22. Liu, G.; Li, C.; Liu, X.; Ge, Q.; Chen, X. Automatic stacking-velocity estimation using similarity-weighted clustering. Geophys. Prospect. 2018, 66, 649–663. [Google Scholar] [CrossRef]
  23. Velis, D. Simulated annealing velocity analysis: Automating the picking process. Geophysics 2021, 86, V119–V130. [Google Scholar] [CrossRef]
  24. Mulder, W.A.; ten Kroode, A.P.E. Automatic velocity analysis by differential semblance optimization. Geophysics 2002, 67, 1184–1191. [Google Scholar] [CrossRef]
  25. Shen, P.; Symes, W.W. Automatic velocity analysis via shot profile migration. Geophysics 2008, 73, VE49–VE59. [Google Scholar] [CrossRef] [Green Version]
  26. Li, J.; Symes, W.W. Interval velocity estimation via NMO-based differential semblance. Geophysics 2007, 72, U75–U88. [Google Scholar] [CrossRef]
  27. Zhang, Y.; Biondi, B. Moveout-based wave-equation migration velocity analysis. Geophysics 2013, 78, U31–U39. [Google Scholar] [CrossRef] [Green Version]
  28. Weibull, W.W.; Arntsen, B. Automatic velocity analysis with reverse-time migration. Geophysics 2013, 78, S179–S192. [Google Scholar] [CrossRef]
  29. Sun, B.; Alkhalifah, T. Automatic wave-equation migration velocity analysis by focusing subsurface virtual sources. Geophysics 2018, 83, U1–U8. [Google Scholar] [CrossRef] [Green Version]
  30. Choi, H.; Byun, J.; Seol, S.J. Automatic velocity analysis using bootstrapped differential semblance and global search methods. Explor. Geophys. 2010, 41, 31–39. [Google Scholar] [CrossRef] [Green Version]
  31. Zuberi, A.; Alkhalifah, T. Imaging by forward propagating the data: Theory and application. Geophys. Prospect. 2013, 61, 248–267. [Google Scholar] [CrossRef]
  32. Zuberi, M.A.H.; Alkhalifah, T. Generalized internal multiple imaging (GIMI) using Feynman-like diagrams. Geophys. J. Int. 2014, 197, 1582–1592. [Google Scholar] [CrossRef]
  33. Singh, S.; Snieder, R.; Behura, J.; van der Neut, J.; Wapenaar, K.; Slob, E. Marchenko imaging: Imaging with primaries, internal multiples, and free-surface multiples. Geophysics 2015, 80, S165–S174. [Google Scholar] [CrossRef] [Green Version]
  34. Leite, L.W.B.; Vieira, W.W.S. Automatic seismic velocity analysis based on nonlinear optimization of the semblance function. J. Appl. Geophys. 2019, 161, 182–192. [Google Scholar] [CrossRef]
  35. Weglein, A.B.; Hsu, S.-Y.; Terenghi, P.; Li, X.; Stolt, R.H. Multiple attenuation: Recent advances and the road ahead. Lead. Edge 2011, 30, 864–875. [Google Scholar] [CrossRef]
  36. Wang, Y.H. Multiple attenuation: Coping with the spatial truncation effect in the Radon transform domain. Geophys. Prospect. 2003, 51, 75–87. [Google Scholar] [CrossRef]
  37. Herrmann, F.J.; Boeniger, U.; Verschuur, D.J. Non-linear primary-multiple separation with directional curvelet frames. Geophys. J. Int. 2007, 170, 781–799. [Google Scholar] [CrossRef] [Green Version]
  38. de Souza, M.S.; Porsani, M.J. Automatic velocity analysis using complex seismic traces. In Proceedings of the 14th International Congress of the Brazilian Geophysical Society & EXPOGEF, Rio de Janeiro, Brazil, 3–6 August 2015; pp. 1337–1341. [Google Scholar] [CrossRef]
  39. He, W.; Geng, J. Elimination of free-surface-related multiples by combining Marchenko scheme and seismic interferometry. Geophysics 2022, 87, Q1–Q14. [Google Scholar] [CrossRef]
  40. Wiggins, J.W. Attenuation of complex water-bottom multiples by wave-equation-based prediction and subtraction. Geophysics 1988, 53, 1527–1539. [Google Scholar] [CrossRef]
  41. Verschuur, D.J.; Berkhout, A.; Wapenaar, C. Adaptive surface-related multiple elimination. Geophysics 1992, 57, 1166–1177. [Google Scholar] [CrossRef]
  42. Verschuur, D.J.; Berkhout, A.J. Removal of internal multiples with the common-focus-point (CFP) approach: Part 2—Application strategies and data examples. Geophysics 2005, 70, V61–V72. [Google Scholar] [CrossRef]
  43. Guitton, A.; Verschuur, D. Adaptive subtraction of multiples using the L1-norm. Geophys. Prospect. 2004, 52, 27–38. [Google Scholar] [CrossRef] [Green Version]
  44. Herrmann, F.J.; Wang, D.; Verschuur, D.J. Adaptive curvelet-domain primary-multiple separation. Geophysics 2008, 73, A17–A21. [Google Scholar] [CrossRef] [Green Version]
  45. Li, Z.-x. Adaptive multiple subtraction based on support vector regression. Geophysics 2020, 85, V57–V69. [Google Scholar] [CrossRef]
  46. Dragoset, B.; Verschuur, E.; Moore, I.; Bisley, R. A perspective on 3D surface-related multiple elimination. Geophysics 2010, 75, A245–A261. [Google Scholar] [CrossRef] [Green Version]
  47. Natan. Fast 2D Peak Finder. 2021. Available online: https://www.mathworks.com/matlabcentral/fileexchange/37388-fast-2d-peak-finder (accessed on 26 May 2021).
  48. Fomel, S. Local seismic attributes. Geophysics 2007, 72, A29–A33. [Google Scholar] [CrossRef]
  49. Hu, B.; Wang, D.; Wang, R. An Iterative Focal Denoising Strategy for Passive Seismic Data. Pure Appl. Geophys. 2020, 177, 4607–4622. [Google Scholar] [CrossRef]
  50. Claerbout, J.F.; Abma, R. Earth Soundings Analysis: Processing versus Inversion; Blackwell Scientific Publications: London, UK, 1992; Volume 6. [Google Scholar]
  51. Fomel, S.; Jin, L. Time-lapse image registration using the local similarity attribute. Geophysics 2009, 74, A7–A11. [Google Scholar] [CrossRef]
  52. Jiménez, A.; Ríos-Insua, S.; Mateos, A. A generic multi-attribute analysis system. Comput. Oper. Res. 2006, 33, 1081–1101. [Google Scholar] [CrossRef] [Green Version]
  53. Ríos-Insua, S.; Gallego, E.; Jiménez, A.; Mateos, A. A multi-attribute decision support system for selecting intervention strategies for radionuclide contaminated freshwater ecosystems. Ecol. Model. 2006, 196, 195–208. [Google Scholar] [CrossRef]
  54. Shevchenko, G.; Ustinovichius, L.; Andruškevičius, A. Multi-attribute analysis of investments risk alternatives in construction. Technol. Econ. Dev. Econ. 2008, 14, 428–443. [Google Scholar] [CrossRef]
  55. Zavadskas, E.K.; Turskis, Z.; Kildiene, S. State of art surveys of overviews on mcdm/madm methods. Technol. Econ. Dev. Econ. 2014, 20, 165–179. [Google Scholar] [CrossRef] [Green Version]
  56. Ferreira, R.S.; Oliveira, D.A.B.; Semin, D.G.; Zaytsev, S. Automatic Velocity Analysis Using a Hybrid Regression Approach With Convolutional Neural Networks. IEEE Trans. Geosci. Remote Sens. 2021, 59, 4464–4470. [Google Scholar] [CrossRef]
  57. Abbad, B.; Ursin, B.; Rappin, D. Automatic nonhyperbolic velocity analysis. Geophysics 2010, 75, Y3. [Google Scholar] [CrossRef]
  58. Van der Baan, M.; Kendall, J.M. Estimating anisotropy parameters and traveltimes in the τ-p domain. Geophysics 2002, 67, 1076–1086. [Google Scholar] [CrossRef]
  59. Sheriff, R.E.; Geldart, L.P. Exploration Seismology; Cambridge University Press: Cambridge, UK, 1995. [Google Scholar]
  60. Chen, C.-T. Extensions of the TOPSIS for group decision-making under fuzzy environment. Fuzzy Sets Syst. 2000, 114, 1–9. [Google Scholar] [CrossRef]
  61. Xu, Z. Uncertain Multi-Attribute Decision Making: Methods and Applications; Springer: Berlin/Heidelberg, Germany, 2015. [Google Scholar] [CrossRef]
  62. Shih, H.-S.; Shyur, H.-J.; Lee, E.S. An extension of TOPSIS for group decision making. Math. Comput. Model. 2007, 45, 801–813. [Google Scholar] [CrossRef]
  63. Zavadskas, E.K.; Mardani, A.; Turskis, Z.; Jusoh, A.; Nor, K.M. Development of TOPSIS method to solve complicated decision-making problems—An overview on developments from 2000 to 2015. Int. J. Inf. Technol. Decis. Mak. 2016, 15, 645–682. [Google Scholar] [CrossRef]
  64. Yoon, K.P.; Kim, W.K. The behavioral TOPSIS. Expert Syst. Appl. 2017, 89, 266–272. [Google Scholar] [CrossRef]
  65. Zhang, Y.; Liu, Y.; Shi, X.; Yang, C.; Wang, W.; Li, Y. Risk evaluation of coal spontaneous combustion on the basis of auto-ignition temperature. Fuel 2018, 233, 68–76. [Google Scholar] [CrossRef]
  66. Vega, A.; Aguaron, J.; Garcia-Alcaraz, J.; Mara Moreno-Jimenez, J. Notes on Dependent Attributes in TOPSIS. In Proceedings of the 2nd International Conference on Information Technology and Quantitative Management (ITQM), the National Research University Higher School of Economics, Moscow, Russia, 3–5 June 2014; pp. 308–317. [Google Scholar] [CrossRef] [Green Version]
  67. Papathanasiou, J.; Ploskas, N. Topsis. In Multiple Criteria Decision Aid; Springer: Berlin/Heidelberg, Germany, 2018; pp. 1–30. [Google Scholar] [CrossRef]
  68. Hua, Z.; Gong, B.; Xu, X. A DS-AHP approach for multi-attribute decision making problem with incomplete information. Expert Syst. Appl. 2008, 34, 2221–2227. [Google Scholar] [CrossRef]
  69. Abdalla, A.; Cen, H.; Abdel-Rahman, E.; Wan, L.; He, Y. Color Calibration of Proximal Sensing RGB Images of Oilseed Rape Canopy via Deep Learning Combined with K-Means Algorithm. Remote Sens. 2019, 11, 3001. [Google Scholar] [CrossRef] [Green Version]
  70. Verm, R.W.; Symes, W.W. Practice and pitfalls in NMO-based differential semblance velocity analysis. In SEG Technical Program Expanded Abstracts 2006; Society of Exploration Geophysicists: Houston, TX, USA, 2006; pp. 2112–2116. [Google Scholar] [CrossRef] [Green Version]
  71. Berkhout, A.J.; Verschuur, D.J. Removal of internal multiples with the common-focus-point (CFP) approach: Part 1—Explanation of the theory. Geophysics 2005, 70, V45–V60. [Google Scholar] [CrossRef] [Green Version]
  72. Ding, C.; Ma, J. Automatic migration velocity analysis via deep learning. Geophysics 2022, 87, U135–U153. [Google Scholar] [CrossRef]
  73. Park, M.J.; Sacchi, M.D. Automatic velocity analysis using convolutional neural network and transfer learning. Geophysics 2020, 85, V33–V43. [Google Scholar] [CrossRef]
  74. Crowe, T.J.; Noble, J.S.; Machimada, J.S. Multi-attribute analysis of ISO 9000 registration using AHP. Int. J. Qual. Reliab. Manag. 1998, 15, 205–222. [Google Scholar] [CrossRef]
Figure 1. Comparison between seismic data without and with multiples and their velocity spectra, There is a corresponding relationship between event and peak marked with the same number in the figure, where the number 1 is primary, number 2 is surface-related multiple, and number is internal multiple: (a) seismic data containing only primaries. (b) Velocity spectra corresponding to data (a). (c) Seismic data containing primaries and multiples. (d) Velocity spectra corresponding to data (c).
Figure 1. Comparison between seismic data without and with multiples and their velocity spectra, There is a corresponding relationship between event and peak marked with the same number in the figure, where the number 1 is primary, number 2 is surface-related multiple, and number is internal multiple: (a) seismic data containing only primaries. (b) Velocity spectra corresponding to data (a). (c) Seismic data containing primaries and multiples. (d) Velocity spectra corresponding to data (c).
Remotesensing 14 05428 g001
Figure 2. A demonstration of the similarity between the original data velocity spectra and the predicted multiple velocity spectra: (a) the velocity spectra of original data. (b) The velocity spectra of predicted multiples. (c) Local similarity spectra.
Figure 2. A demonstration of the similarity between the original data velocity spectra and the predicted multiple velocity spectra: (a) the velocity spectra of original data. (b) The velocity spectra of predicted multiples. (c) Local similarity spectra.
Remotesensing 14 05428 g002
Figure 3. Local similarity of the picked peaks.
Figure 3. Local similarity of the picked peaks.
Remotesensing 14 05428 g003
Figure 4. Display of synthetic data: (a)Velocity model containing a velocity reversal, (b) CMP record of original data, (c) velocity spectra of original data.
Figure 4. Display of synthetic data: (a)Velocity model containing a velocity reversal, (b) CMP record of original data, (c) velocity spectra of original data.
Remotesensing 14 05428 g004
Figure 5. Display of synthetic data test results: (a) muting-based strategy (MU), (b) Radon transform strategy (RT), (c) multiple attenuation strategy (MA), (d) proposed advanced method (AD).
Figure 5. Display of synthetic data test results: (a) muting-based strategy (MU), (b) Radon transform strategy (RT), (c) multiple attenuation strategy (MA), (d) proposed advanced method (AD).
Remotesensing 14 05428 g005
Figure 6. Velocity comparison.
Figure 6. Velocity comparison.
Remotesensing 14 05428 g006
Figure 7. NMO comparison results. 1.4–2.4 s: (a) NMO results of the MP velocity, (b) NMO results of the AD velocity, (c) NMO results of the MA velocity. 3.2–4.2 s: (d) NMO results of the MP velocity, (e) NMO results of the AD velocity, (f) NMO results of the MA velocity.
Figure 7. NMO comparison results. 1.4–2.4 s: (a) NMO results of the MP velocity, (b) NMO results of the AD velocity, (c) NMO results of the MA velocity. 3.2–4.2 s: (d) NMO results of the MP velocity, (e) NMO results of the AD velocity, (f) NMO results of the MA velocity.
Remotesensing 14 05428 g007
Figure 8. Application on a single CMP gather from the field data: (a) CMP record of original data, (b) velocity spectra of original data, (c) velocity spectra and velocity curve of the control group, in which the red solid line is automatically obtained, and the yellow dotted line is the reference velocity, (d) velocity spectra and velocity curve obtained by the method proposed in this paper, in which the red solid line is automatically obtained, and the yellow dotted line is the reference velocity.
Figure 8. Application on a single CMP gather from the field data: (a) CMP record of original data, (b) velocity spectra of original data, (c) velocity spectra and velocity curve of the control group, in which the red solid line is automatically obtained, and the yellow dotted line is the reference velocity, (d) velocity spectra and velocity curve obtained by the method proposed in this paper, in which the red solid line is automatically obtained, and the yellow dotted line is the reference velocity.
Remotesensing 14 05428 g008
Figure 9. Comparison of velocity model: (a) reference model, (b) velocity model obtained by control group, (c) velocity model obtained by proposed method.
Figure 9. Comparison of velocity model: (a) reference model, (b) velocity model obtained by control group, (c) velocity model obtained by proposed method.
Remotesensing 14 05428 g009
Figure 10. Comparison of stacked sections: (a) stack result of reference velocity, (b) stack result of the velocity obtained by control group, (c) stack result of the velocity obtained by proposed method.
Figure 10. Comparison of stacked sections: (a) stack result of reference velocity, (b) stack result of the velocity obtained by control group, (c) stack result of the velocity obtained by proposed method.
Remotesensing 14 05428 g010
Table 1. Local similarity statistics between primaries and multiples.
Table 1. Local similarity statistics between primaries and multiples.
MaximumMinimumAverageMedian
Primary6.0 × 10−61.6 × 10−112.1 × 10−63.0 × 10−7
Multiple2.1 × 10−14.2 × 10−39.7 × 10−27.8 × 10−2
Table 2. The difference between the reference and estimated velocities.
Table 2. The difference between the reference and estimated velocities.
MBRTMSMI
σ 6.0%5.4%5.0%1.3%
Table 3. Analysis for the color spectra.
Table 3. Analysis for the color spectra.
TypeFeaturesColor
Primarylow similarity, high velocity and amplitude in commonRemotesensing 14 05428 i001
Internal multiplelow amplitude, high similarity. velocity relates to the developed formation.Remotesensing 14 05428 i002
Surface-related multiplehigh similarity, velocity relates to shallow formation, amplitude could be high.Remotesensing 14 05428 i003
Suspicious peakthe peak should be focused onRemotesensing 14 05428 i004
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Zhang, J.; Wang, D.; Hu, B.; Gong, X. An Automatic Velocity Analysis Method for Seismic Data-Containing Multiples. Remote Sens. 2022, 14, 5428. https://doi.org/10.3390/rs14215428

AMA Style

Zhang J, Wang D, Hu B, Gong X. An Automatic Velocity Analysis Method for Seismic Data-Containing Multiples. Remote Sensing. 2022; 14(21):5428. https://doi.org/10.3390/rs14215428

Chicago/Turabian Style

Zhang, Junming, Deli Wang, Bin Hu, and Xiangbo Gong. 2022. "An Automatic Velocity Analysis Method for Seismic Data-Containing Multiples" Remote Sensing 14, no. 21: 5428. https://doi.org/10.3390/rs14215428

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop