Inversion of ocean-bottom seismometer (OBS) waveforms for oceanic crust structure: a synthetic study

The waveform inversion method is applied—using synthetic ocean-bottom seismometer (OBS) data—to study oceanic crust structure. A niching genetic algorithm (NGA) is used to implement the inversion for the thickness and P-wave velocity of each layer, and to update the model by minimizing the objective function, which consists of the misfit and cross-correlation of observed and synthetic waveforms. The influence of specific NGA method parameters is discussed, and suitable values are presented. The NGA method works well for various observation systems, such as those with irregular and sparse distribution of receivers as well as single receiver systems. A strategy is proposed to accelerate the convergence rate by a factor of five with no increase in computational complexity; this is achieved using a first inversion with several generations to impose a restriction on the preset range of each parameter and then conducting a second inversion with the new range. Despite the successes of this method, its usage is limited. A shallow water layer is not favored because the direct wave in water will suppress the useful reflection signals from the crust. A more precise calculation of the air-gun source signal should be considered in order to better simulate waveforms generated in realistic situations; further studies are required to investigate this issue.


Introduction
Seismic waveform inversion is a technique to extract quantitative information on subsurface structure by fitting the synthetic seismograms with that of observation. Through comparison with a travel-time-based method, fullwaveform inversion can improve the resolution of the model as it uses both the amplitude and phase information of various seismic phases contained in seismograms. As the wave equation is directly solved in this method, higherorder effects such as diffractions and multiple scattering are accounted for automatically (Pratt 1999). Many studies on the regional continental crust and upper-mantle structure have been carried out using waveform inversion (Das and Nolet 1998;Sherrington et al. 2004;Shibutani et al. 1996;Zheng et al. 2006;Zhu et al. 2006). Full-waveform inversion is a strong nonlinear optimization process. For gradient-based local searching methods, the final model is strongly related to the initial model and may be trapped into local minima. It is computationally time-consuming to perform large-scale forward modeling of seismic wave field propagation (Virieux and Operto 2009).
For regions where subsurface structure can be approximated by a layered model, global optimization methods, such as genetic algorithm and simulated annealing, can be applied to full-waveform inversion to improve its efficiency and convergence (Sen and Stoffa 2013). A layered 1D velocity structure is the most simplified, but fundamental, approximation of the Earth. It plays a significant role in accurate seismic location, calculation of Green's function, and can serve as the initial model for inversion of fine 2D or 3D velocity structure. A more accurate layered crust model could also contribute to improved study of the deeper Earth structure. Recently, waveform inversion has been conducted at several regions to extract layered continental crustal and upper-mantle structure using the data recorded from local seismic networks with global optimization methods (Abdelwahed and Zhao 2014;Chang and Baag 2006;Li and Lei 2014a, b;Li et al. 2007Li et al. , 2012. The oceanic crustal and upper-mantle structure is less well investigated than continental cases because of lack of seismic observation data. Researchers attempted to invert waveform data obtained from ocean reflection experiment for elastic parameters of oceanic crust (Igel et al. 1996;Mendes et al. 1990) and crust and upper-mantle velocity structure from surface waves (Cara and Lévêque 1987;Debayle and Lévêque 1997). In recent years, ocean-bottom seismometer (OBS) techniques have provided increasing amount of seismic data to study the oceanic crust of localized areas, such as north-eastern Japan (Takahashi et al. 2004), Lucky Strike segment (Seher et al. 2010), Southwest Indian Ridge (Li et al. 2015) and north-eastern South China Sea (Zhao et al. 2010). Full-waveform inversion has been applied to OBS data to improve the resolution of oceanic crust structure (Jian et al. 2014;Operto et al. 2006;Borisov and Singh 2015).
This paper attempts to extend a full-waveform inversion scheme used for continental crust and upper-mantle structure studies to the investigation of the oceanic crust using OBS data. Synthetic OBS data generated from an air-gun source in water is used. One significant difference between land and OBS observation is the presence of the water layer, which produces multiples and strong direct waves that cannot be used in waveform inversion. The validity and efficiency of OBS data methods are studied with the use of synthetic numerical models. Considering a limited number and sparse distribution of OBS stations, compared with the land seismic network, an attempt is made to perform inversion for a layered oceanic crust model. The parameters of interest include thickness and P-wave velocity of the sedimentary layer and crust. A niching genetic algorithm (NGA) is employed to perform the global optimization for the waveform inversion. The influence and performance under different observation systems of some key parameters in the NGA method are investigated, and a strategy is proposed to accelerate the convergence rate without increasing the computational complexity. This strategy is verified through comparative tests.

Waveform inversion method
Waveform inversion is a highly nonlinear problem. There are many traditional methods that can be used to determine seismic velocity structure through the use of waveform inversion, such as the conjugate gradient method, the grid search method, genetic algorithms (GA), and simulated annealing. Their common point is that they can only give one minimum of the objective function as the final model, and it is possible that the solution may converge to a local minimum if the initial model is quite far away from the real model (Maurice et al. 2003). However, NGA can search different minima by simulating the evolutionary processes in biology, such as crossover, mutation, selection, and competition (Mahfoud 1995). Hence, NGA is quite suitable for a multimodal optimization problem in geophysics. Koper et al. (1999) developed the NGA and applied it in a teleseismic waveform inversion for the source parameters of the M w 7.2 Kuril Islands earthquake in 1996. Maurice et al. (2003) inferred the crustal and upper-mantle structure under southernmost South America through the use of NGA. Lawrence and Shearer (2006) gave a constrained seismic velocity and density for the mantle transition zone, and Li et al. (2012) inferred the 1D crustal structure under south-eastern Gansu, China, by applying NGA to regional waveform inversion.
Several researches have investigated the effects of different objective functions and parameters in NGA, such as the number of models in each subpopulation, and the critical separation radius (Koper et al. 1999;Li and Lei 2014a). In this paper, further discussion of NGA implementation is presented based on these previous studies, and an attempt is made to illustrate some other factors that could influence the resultant final model.
The NGA works as follows: first of all, initial models are created and then divided into n subpopulations, or demes. A deme represents a group of models that are close to each other in the whole model space, and there is a distance between different demes. We set each deme to contain m models. Traditional genetic algorithm approach is carried out in each deme to work out the second deme, but in NGA the similarity of each member of this deme is calculated with respect to the best model from the former ones. Similarity is an important parameter. The simplest approach is to define the distance, D, of two models, x and y, as the arithmetic average of the normalized separation of the model parameters (Koper et al. 1999), as shown in Eq. (1), where m represents the number of parameters we need to search, x i and y i are the ith parameters of models x and y, and b i and a i are the upper and lower bounds of the ith parameters. The similarity varies from 0, for two identical models, to 1 for two models at opposite ends of the search boundary. If one's similarity exceeds a specified criterion (called R c ), then that model is given a high penalty and is eliminated in the next generation. This process continues until it reaches the final generation.

Model parameterization
The waveform inversion assumes a layered, laterally homogeneous model and seeks to invert for crustal and upper-mantle thickness and P-wave velocity. The S-wave velocity for each layer is derived from the relationship v S = kv P , where k can be determined by the Poisson ratio of each layer. For the crustal layer, k is set to 0.577, which means that the Poisson ratio is equal to 0.25 for this case.
For the sediment layer, k is equal to 0.489. The synthetic model is designed by considering the oceanic crust model given by Rao et al. (2012) and Arnulf et al. (2014). The model parameters are listed in Table 1.

Objective function and synthetic seismograms
The primary component of the objective function is the average value of root-mean-square residual between the synthetic and observed seismograms and their cross-correlation in the time domain. This is shown in Eq. (2): where OðtÞ Â SðtÞ ¼ R OðsÞSðt À sÞds, and N w is the number of waveforms used. The first term of this formula is the root-mean-square error, O ij and S ij are the amplitudes of the observed and synthetic seismograms of the ith component at the jth sampling point. The second term contains the cross-correlation of the synthetic and observed waveform, which shows their similarity. O i and S i are the ith waveform of the observed and synthetic seismograms. This objective function takes into account both the amplitude and phase information of each seismic phase (Li et al. 2012).
A seismic wave field generated from an air-gun source in the water is studied. Synthetic waveforms are calculated using a reflectivity method (Fuchs and Müller 1971). In comparison with land-based observation, a major difference to OBS data is the existence of the water layer, which produces multiples in the water and strong direct waves that cannot be used in the inversion for subsurface structure. In the forward modeling, the water multiples are removed by omitting them from the calculation since they are useless for the crustal structure, though they are of a high amplitude. The direct wave propagating in the water cannot be ignored, and poses the biggest challenge to this study. The very strong direct wave and relatively weak series of reflection waves from the crust are visible in Fig. 1. With the increase of the offset, the direct wave tends to become mixed up with the first reflection phase. This prevents the study of a rather large area, restricting the study region in which the direct wave and reflection wave are separated. The useful reflection signals for inversion are selected, and the direct wave is ignored.  The seismic source is an explosive type, and a Ricker wavelet is used to imitate the air-gun signal. The source is located at a depth of 20 m in the water layer. Considering that the real air-gun signal contains more high-frequency components, a Butterworth band-pass filter is applied to the source with corner frequencies of 5 and 15 Hz, and a central frequency of 10 Hz.
The first phase with a large amplitude shown in Fig. 1 is the direct wave, which adds a large noise signal to the reflected waves. Since they are not overlapped, a time window can be used to pick up the subsequent reflections for waveform inversion.

General test
As shown in Table 1, a total of five P-wave velocities and four thicknesses of layers are searched. The searching ranges are within ±0.5 km/s for v P and ±30 % for the thickness around the preset values described in Table 1. In the first test, ten receivers are set, in line, at the surface of the sediment layer, each between 100 meters and 5 km in horizontal distance apart from the source, as shown in Fig. 2. The NGA inversion is calculated for 500 generations with ten demes, with each deme containing ten models. Seven different NGA runs are performed, using different random model generators (by changing the random number value in the NGA). The average of seven results is taken, with the minimum value of the objective function in each run being used for the final model, with their standard deviation as the error. The probabilities of crossover and mutation are set to 0.90 and 0.10, respectively.
The results shown in Table 2 and Fig. 3 show the comparison of the inversion result with the true model. Figure 4 shows the cost-generation curve. It is clear that the convergence rate slows down gradually, and that the change becomes negligible after hundreds of generations. It is evident from Fig. 5 that the waveform is a good fit.

Number of models in each deme
The number of models (N m ) in each deme has a great effect on the result of the NGA inversion. If N m is too small, then the rate of convergence will be low if there is poor diversity in each subpopulation. This could easily be solved by increasing the number of models, which would also lead to more computing time. So there is a trade-off between improving the convergence rate and cutting down the computing cost. A test is performed with N m equal to 6, 8, 10, 12, and 16. It is evident from Fig. 6 that the results of the test with the smallest value of N m have the lowest quality, and that quality improves as N m increases. It was found that once N m has reached a value greater than 10, then there is no further apparent change in the cost value. For an appropriate balance of accuracy and time requirement, N m equal to 10 was selected.

Critical distance R c
It is always the hope that the inversion method can search the whole solution space and avoid converging on the local minima as much as possible. In NGA, this is carried out by choosing a suitable R c value, which controls how different subpopulations migrate into different niches of the solution space. If R c equals 0, then all demes will inhabit global minimum as no artificial distinction exists between each deme, like performing several GA at the same time. However, if R c is too large, then only the best model evolves with time, while the others are artificially ruled out and become ''frustrate''. According to Koper et al. (1999), the cost value of a model that is ''frustrate'' always remains high and shows little improvement over time, being reinitialized randomly after every generation. This is undesirable behavior, because competition and selection between useful demes is required in order to efficiently search for the optimum model.
Four different R c values are selected for this test: 0.01, 0.1, 0.2, and 0.3. Values greater than 0.3 resulted in poor performance in pretest and are not discussed. It is evident that the cost value is not so favorable when R c is equal to 0.3, and from Fig. 7 there is no obvious distinction between the other tests. It is only by examining the results of the 200th generation that some subpopulations are separated when R c equals 0.1 or 0.2. With R c equal to 0.01, most subpopulations converge to the global minimum. Considering the competition between subpopulations, it is desirable to set R c between 0.1 and 0.2 in order to make our NGA different from the traditional GA, while at the same time ensuring that none, or very few, subpopulations become ''frustrate''-ensuring the search efficiency of the NGA.

Effect of different observation systems
In the previous discussion, the receivers were always arranged in a line, with the goal of achieving an along-axis profile. However, there are several ways to set the observation system based on different research purposes. This section discusses the feasibility of NGA inversion methods using a variety of observation systems.    Fig. 8. In some cases, it was not possible to set a sufficient number of OBS in line, due to the constraints of the situation. Furthermore, data from the irregular distribution of OBS give the average 1D structure of a wider area, while linear arrangement just focuses on the profile. Therefore, it is necessary to perform tests for such an OBS network. The results of these tests are presented in Table 3 and Figs. 9, 10, and 11. A good fitting of the waveform is still observed.

Single receiver
Even though the observation system could include many accurately timed, high-gain seismometers, a single-station method still has advantages, of which economic cost is an important one. Tests were performed in order to examine whether single OBS data are sufficient to explore the local crust structure. The distance of the OBS to the source is 3.5 km. The results of the tests are presented in Table 4 and Figs. 12, 13, and 14. The waveform fitting is sufficiently good to validate the single receiver method.

Strategy to speed up convergence
Through observation of the evolution of cost values in the previous numerical tests, it is apparent that calculating for more generations gives a final model with a lower cost value, but that this comes at the expense of significantly increased time consumption. Because significant changes rarely take place after about 100 generations, as shown by the cost-generation curves, it is proposed to use the result of the first few generations to further narrow the range of each parameter to conduct a second inversion. Generally speaking, a smaller search range leads to a better inversion result, and this can be judged by comparison of the value of the object function. In this section, the two experiments described in Sect. 5 are revisited, with the results of the 15th generation used. The mean value and standard deviation r of each parameter are calculated, and used to set the upper and lower bound of a new search range equal to mean ±3r. A coefficient smaller than 3 would be more favorable, but would increase the risk that the parameter value of the real model would not be included in the search range, thus making this inversion a failure. For this reason, 3 is chosen as the coefficient, which gives a limited restriction on the range as shown in Tables 5 and 6.  Tables 5 and 6, it can be seen that the range restriction greatly limits the range of the top two layers, especially the first layer, and that improvements on the upper layers have significant consequences. Because each reflection wave phase travels through upper layers, if the velocity and thickness values deviate away from the actual ones too much, then it will be difficult for all of the seismic phases to arrive at the correct time, which increases the difficulty of waveform fitting. Therefore, better estimation of the upper layers is crucial, and the range restriction strategy contributes to a reduction of its uncertainty.
The essence of this strategy is to reduce the searching range of model parameters from the first inversion within a Fig. 9 Comparison of the average model (red line) from seven NGA inversions with the true models (blue line). The preset ranges for velocity and depth in which the model parameters are allowed to change in the inversion are shown by black lines, for the tests described in Sect. 5.1 Fig. 10 Cost versus generations. This result is from the situation in which OBS are located irregularly and sparsely, for the tests described in Sect. 5.1  few generations, and then to carry out a new, second, inversion with the narrowed range. It is seen that convergence on the global minimum is more likely after several generations in the second inversion. By applying this strategy, it takes less than 50 generations to obtain a result that has a similar cost value as was found from the ordinary inversion as shown in Fig. 15. A time saving of almost 80 % is made in this way, demonstrating that the computational efficiency is greatly enhanced. This strategy utilizes a simple ''mean ± kr'' approach, but there remains a question of optimal selection of the constant k. This work set k equal to 3, and the results show Fig. 12 Comparison of the average model (red line) from seven NGA inversions with the true models (blue line). The preset ranges for velocity and depth, in which the model parameters are allowed to change in the inversion, are shown by black lines, for the tests described in Sect. 5.2 Fig. 13 Cost versus generations. This result is from the situation in which the OBS is located 3.5 km away, for the tests described in Sect. 5.2  that 3 is a suitable value, but it is not proven to be the ideal value. It may be difficult to prove the existence of an optimal k value, and so the value should be based on practical experience of real cases.

Discussion and conclusions
This work attempted to invert the waveform recorded by OBS using a synthetic model, for the purpose of studying the structure of the oceanic crust. Tests were carried out on the parameters of a niching genetic algorithm in order to determine appropriate values. Parameter values of N m equal to 10 and of R c between 0.1 and 0.2 were found to be favorable. Different observation systems, irregular distribution and single station, were studied, and results showed that both offered quite good waveform fitting. A strategy was proposed to accelerate the convergence rate and its effectiveness was verified. In conclusion, using OBS data for waveform inversion is effective in the study of a Fig. 15 Comparison of cost-generation curves of prerange-restriction test (red) with those with the range restriction strategy applied (blue): a the irregular and sparse OBS distribution described in Sect. 5.1 and b the single OBS test described in Sect. 5.2. The sharp rise is caused by the reshuffling of all the parameter values from the new, restricted ranges to determine the initial model for the second inversion 5.00-6.00 5.02-5.96 6.7 6.20-7.20 6.23-7.14 7.1 6.60-7.60 6.60-7.60 8.1 7.60-8.60 7.60-8.60 Earthq Sci (2016) 29 (4):203-213 211 regional oceanic crust structure that can be approximated by a layered model. Limitations as well as advantages were found. Most significantly, this method is confined to a circular range of several kilometers, where the direct wave in water does not overlap with the useful reflection signals. By setting the thickness of the water layer to 3.2 km, a circular area of 5 km radius is available. This radius will be smaller if the depth of the sea level is shallower, as the direct wave will suppress the reflected wave at a close distance. Measures should be taken in order to study larger areas, and in particular, an effective method is required to ignore the direct wave while not adding artificial noise to the reflection signals. In addition, this study used a Ricker wavelet to approximate the air-gun source signal, but the real seismic source is much more complicated. Many researchers have discussed improved approximation of the real bubble behavior (e.g., Johnson 1994;Schulze-Gattermann 1972), but this was not included in this study. This study did not consider the influence of different source time functions, but if an accurate expression of air-gun source time function is calculated, then this will offer a better simulation of the real situation.