{'fraction_non_alphanumeric': 0.02746427333307654, 'fraction_numerical': 0.0077038634875390005, 'mean_word_length': 4.8028609745194455, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Comparing performance of multi-frequency bands Occam’s receiver function inversion to standard linearized receiver function inversion

Receiver function (RF) is a characterized waveform that sensitive to the sub-surface velocity structure near the seismic station. In the geophysical problem, an inversion scheme has been used for determining the structure model that corresponded to the observed data. Because RF inversion is a non-linear problem, the standard method such as linearized inversion suffered in the non-uniqueness of the results. Thus, using a stochastic optimization algorithm (e.g. Monte Carlo or Generic Algorithm) is the most interesting trend for solving the problem. However, their computation cost is high and the efficiency is strongly depended on user setting. In this work, Occam’s inversion algorithm, which is popular in the electromagnetic method, and multi-frequency bands receiver function have been integrated to improve the efficiency of receiver function inversion. The propagation method has been applied for the calculation of seismogram in this work. The synthetic cases, which are thin low velocity zone model (TLVZ), broad low velocity zone model (BLVZ), thin high velocity zone model (THVZ) and broad high velocity zone model (BHVZ), were used to compare the performance of this newly implemented algorithm and the standard linearized inversion. Single frequency band (SFO) and multi-frequency bands Occam’s inversion (MFO) have been tested for each case. For, TLVZ, THVZ and BHVZ, the MFO inversion provided better fitting results than SFO in most cases. Comparing to linearized inversion of CPS330, the MFO results still provide better performance. For the smooth low velocity structure in BLVZ model, both linearized and MFO can recover a significant structure of the true model. More concisely, the linearized inversion can recover the absolute velocity of the true model better than MFO. By the way, the MFO provides the best data fitting for the BLVZ. In summary, the MFO receiver function inversion that implemented in this study provides a new improved tool for a seismologist to investigating the subsurface structure.


Introduction
The seismic wave generated by earthquake contains information of source characteristic, earth structure along the propagation path, and the response of the instrument. The strong signal of seismic waves from the earthquake gives an advantage for studying the earth's interior. Among many seismological  2 methods, one of the novel seismological methods is the receiver function (RF) method [1]. The plane wave assumption of teleseismic P wave provided a useful assumption that can be used to extract the information of the earth structure. For the receiver function method, the conversion and reverberation phases of P-wave can be used to characterize the crustal structure beneath a seismic station. For example, a two layers crustal model and its RF response are presented in figure 1. RF is more complicated for real structure because there are many other discontinuities within the earth such as Conrad discontinuity at mid-crust, sedimentary, magmatic underplating, or mantle transition zone. The complex earth structure makes receiver function time series complicate and difficult to be interpreted visually.
Naturally, the earth velocity structure control the characteristic of the observed receiver function (figure 1). Vice versa, in geophysical research, its inverse problem is to mathematically convert an observed receiver function time series into a velocity model of the unknown subsurface structure. Because of the non-linear characteristic of the problem, the inverse problem of the receiver function is very challenging. The linearized receiver function inversion [2] has shown that receiver function inversion can be suffered from the non-uniqueness of results. Several stochastic optimization algorithms have been introduced to solve receiver function inversion (e.g. Malinverno, 2010) [3]. However, those algorithms are largely depending on user-defined parameters and high computational cost. In this study, two approaches will be implemented into receiver function inversion. First approach is Occam's algorithm [4]. It has been effectively used in geophysical problem especially in the geo-electromagnetic method [5]. The core idea of Occam's method is to make the smoothest model that can produce a good fitting response to the data. The second approach is the multi-frequency bands receiver function inversion. Bodin (2012) [6] and Li (2017) [7] shown the improving performance of inversion by applied multi-frequency bands data to the receiver function stochastic optimization. High-frequency response in RF is important for the detection of sharp velocity structure, and the low-frequency response contains information about the small wave number structure or the gradually change velocity structure [8]. The incorporation of multiple-frequency bands RFs will increase the accuracy of an inverted model from the inversion.

Methodology
This paper focuses on the performance testing of newly implemented receiver function inversion technique, the multi-frequency bands Occam's receiver function inversion. This inversion combined the advantage of Occam's inversion and multi-frequency bands receiver function to enhance the solvability of the prospective velocity model. The first two subsections briefly describe the important parts of the inversion algorithm. The last subsection describes the strategies and criteria that use to analyze inversion performance.

The Occam's inversion
The meaning of inverse problem is to find a suitable model that matches with the objectives of the problem. Generally, in geophysical problem, the main objective is to construct a model that produces a predicted or calculated data that fits well with the observed data. To avoid an unrealistic solution, the second objective of model constraint was introduced to govern the range of solution in the inverse problem. The function that combined these two objectives is then called "objective function". In this work, an objective function (U) follows the modified version of Occam's equation from equation 3.5 of Vachiratienchai [9]. . (1) Here, the first term on the righthand side of equation (1) is the model constraint. It is the difference between the model (m) and the prior model (m 0 ) that characterized by the model covariance matrix, C m . The second term presents the misfit term, which is the difference between observed data (d) and predicted data (F[m]). The data covariance matrix (C d ) obtained from the diagonalized matrix of reciprocal of standard error of data. The Lagrange multiplier (λ) is a parameter that controls the tradeoff between model constraint and data misfit. For the calculation of predicted data from the model, the 3 receiver function forward operator (F) included two steps. In the first step, radial (R) and vertical (V) seismogram are calculated by using the propagation matrix method [10]. Second, to obtain the predicted receiver function, the radial seismogram ( ) is deconvoluted by vertical seismogram ( ) in the frequency domain using a water-level method [1]: where is a complex conjugate of and G is introduced as a low-pass Gaussian filter. The parameter α is the Gaussian width parameter. The is water level parameter that used to avoid the zero value of the denominator in the frequency domain. Where c is the lowest floor of the denominator. The suitable water level value that is the smallest value that produced smallest artificial noise and stable division result.
The goal of the receiver function inverse problem is to find the model that minimized the objective function in equation (1). The optimization form of objective function after differentiating respect to the model parameter (m) can be written in the iterative form [9] as equation (5).   (5) where and . J is Jacobian matrix or sensitivity, which represent the changing in predicted receiver function respect to changing in the model, and the subscription denoted the iteration number. In the Occam's algorithm, the Lagrange multiplier will be searched for the most suitable value by regularization method [5]. To obtain the best fitting model, an inversion algorithm (equation 5) is performed iteratively until the objective function is minimized or the model is converged. The flow chart of an inversion procedure is shown in figure 2.

Figure 2. Occam's inversion procedure flow chart
For the first iteration, the starting prior model is set to be the same as the initial model, which is defined by user assumption. Although, the initial model can be derived from the geological or other geophysical constraints. In each iteration, when the new inverted model is calculated, but the new misfit is not decreased from the previous iteration, the prior model will be replaced by latest inverted model. Otherwise, the current inverted model will be the initial model for the next iteration. These processes will run iteratively until two criteria were met. First, the misfit is smaller than the target misfit, and the second, the iteration number is equal to maximum iteration. Both target misfit and the maximum iteration number are defined by the user.

Parameterization and multi-frequency bands approach
In this inversion, the model parameter (m) is the one-dimension horizontal layers of S-wave velocity (V s ). While the thickness and V p /V s ratio of each layer are fixed to reducing the computational cost. The P-wave velocity model is calculated by V p /V s ratio. The thickness used in section 2.3 is 1.25 km same 5 as a reference model in Ammon [2]. The density model (ρ) given by Birch's relation [11] is The water level for the deconvolution process was set as minimum 0.0001.
Because of the difficulty in the noise level measurement of the receiver function time series, the standard error (σ) has been set as 3% of the maximum amplitude. To enhance the result by using multi-frequency bands approach, the observed data (d) is the combination of filtered receiver function from different frequency bands. Hence, the observed multifrequency bands receiver function data can be written as: (6) when d fi is vector of the filtered receiver function at different gaussian factor (number i th ) in equation (3). The parameter n is the number of frequencies data using in the inversion process.

Synthetic test and performance comparison
In this study, four velocity models were used to generate synthetic data that assumed to be an observed receiver function. The velocity models using in the test follow the synthetic test of Ammon [2] who showed the problems on the non-uniqueness of receiver function inversion. Here, all of four models include the structure from surface to a depth of 50 km. The model is divided into multiple layers with 1.25 km thick. The S-wave velocity in the crust is starting from 3.4 km/s at the surface. Then, the velocity is increased with depth until reach a sharp velocity increase at the mantle, the so-called "Mohorovicic discontinuity or Moho". The velocity of the mantle in these four models is 4.6 km/s. For the first model, thin low-velocity zone (TLVZ), the sharp 2.5 km thick of low velocity layers (3.2 km/s) were added to the crust at depth of 20 km. The Moho depth was set to be 32.5 km. The second model, broad-low velocity zone (BLVZ), the slowly decrease low velocity layers were added to the middle crust. While its low velocity zone has a higher velocity than the TLVZ model, it contains a broader zone of the low velocity layers. The BLVZ model has a thicker crust and Moho depth is 35 km. The third model, thin high velocity zone (THVZ), is comparable to the TLVZ model, but a thin low velocity zone was replaced by the thin high velocity zone. The fourth model, broad high velocity zone (BHVZ), opposed to the BLVZ model where a low velocity zone replaced by a high velocity zone and Moho depth of this model is equal to first and third model. These four models are shown as "true" models in figure 3. The synthetic seismogram generated from these models will be used as observed data for inversion. In this test, no noise was adding to the data. Four synthetic data set from models mentioned above were applied to the inversion algorithms. For each model test, three different inverted models from three inversion methods will be compared. First, the single frequency band Occam's receiver function inversion (SFO) will be a candidate to present the performance of Occam's inversion. The frequency using for SFO is receiver function with Gaussian factor equal to 2.5 (cut-off frequency ~ 1Hz). This frequency is high enough to cover the smallest anomaly using in the synthetic test. Second, the multi-frequency bands Occam's inversion (MFO), four filtered receiver function waveforms from different frequency bandpass were used and combined as equation (6). The Gaussian widths used in this study are 1.0, 1.5, 2.5, and 4 which can recover both sharp and broad structure. The last, hrftn96 program in Computer Program in Seismology, CPS330 [12] is used as a competitor and reference algorithm for the inversion. The hrftn96 is based on the linearized inversion of Ammon [2] and has been popularly used by many seismologists. Similar to MFO, the multifrequency bands data will be applied to CPS330 inversion. The result of CPS330 and MFO then will present the difference between the Occam's inversion and the linearized inversion. The initial models applied to four model tests and three inversion algorithms are similar. The crustal S-wave velocity is 3.82 km/s. The mantle S-wave velocity is 4.65 km/s. The Moho depth is 32.5 km, and V p /V s ratio is 1.7.
The performance of three RF inversion approaches will be compared using both qualitative and quantitative parameters. where N m is number of model layers. The value of ϕ M quantitatively presents the reliability of results. The lower value of model difference means the results is closer to the correct solution and the algorithm has better performance. However, in practice, the model difference cannot be calculated because the true velocity model of the earth is unknown. The second quantity for the performance analysis is misfit parameter which expresses the difference between predicted data from inversion and observed (or synthetic) data. The misfit is common for inverse problem for the justification of result. In this study, the misfit has been calculated as a root mean square difference between synthetic data (d) and predicted data F[m] as: where N d denotes the number of the data and the index j denoted the data index. The value of misfit presents the goodness of inversion. It quantitatively shows the similarity between the predicted data and the data. The more similarity leads to the lower misfit value. Although, due to the non-uniqueness of the problem, a poor inverted model is also able to produce a low misfit value. Normally, these two quantitative parameters are enough for the comparison of inversion results. Although, there is an ambiguity with in these two values. In the similar way with Ammon [2], The human based interpretation is also needed for the performance analysis. In this study, three other qualitative indexes are applied to the inversion results for the simplification of the comparison. The performance ranking of qualitative indexes is evaluated as A, B and C from best performance, respectively. The detail of evaluation is described as followed.
First, the Moho recovery rank is used for the interpretation of result related to the Moho feature. The Moho is a significant feature that provides the dominant signal in RF. It is also an important constraint for the study of lithospheric structure. Generally, Moho is the easiest feature to be recovered in the receiver function inversion. However, it can be distorted by the complexity of other layering structure. For example, the negative reverberation peaks generated by Conrad discontinuities can be superposition with the main conversion peak from Moho. The inversion algorithm is then unable to solve for the correct Moho features. To consider this Moho significant, the Moho recovery index is using for the comparison of depth and sharpness of the velocity contrast at Moho. First, we consider the depth of the inverted model's Moho and true model's Moho. The method that can provide a correct or a closer value of Moho depth will have a greater rank. If those methods provide the same Moho depth, the velocity contrast of Moho will be compared for the further ranking of this parameter.
Second, the absolute velocity in the inverted models is considered. Ammon [2] shown that, in many cases, inversion can recover only the shape of the anomaly-feature but not the true velocity values. The inverted velocity model is then composed of two velocity characteristics, the shape of the velocity model (relative velocity) and the true velocity or the average trend of the velocity model (absolute velocity). For example, the model with a good relative velocity structure but poor absolute velocity values can provide a small misfit value which might lead to the misinterpretation. Because of the different velocity in crust and mantle, the absolute velocity in the inverted model is considered from the trend of average velocity in crust and mantle, separately.
Third, the shape of the predicted waveform is another characteristic to be considered. The misfit value provides the first order information on the similarity between predicted receiver function and synthetic receiver function. However, its value can be dominantly affected by only the main phase from strongest velocity contrast such as Moho. The consistency of small phases in the predicted receiver function and data are then ignored by the misfit values. To assess the existence of small phases, the predicted data is visually compared with the data by the counting of the aligned phases between the predicted data and synthetic data.

Results and discussions
For each synthetic data set from each model in section 2.3, the inverted models, true models, and initial models were plotted together in each panel of figure 3. Figure 4 shows the synthetic receiver function from the true model and the predicted receiver functions from the inverted models from each method. A summary of inversion performance is shown in table 1. Table 1. Comparisons of the inversion performance. Misfit value is fitting between synthetic data and predicted receiver function. The Moho indication means a performance of inversion to recover velocity contrast and depth of Moho. Anomaly means the anomaly recovering for both velocity and depth. The waveform shape is the overall trend of predicted receiver function compare to synthetic receiver function. This can be implied by the misfit value. The model difference shows the summation of velocity difference between the true model and the inverted model. The letters A, B, and C represent the performance quality from the best performance, respectively.

Model
Performance index

Thin low velocity anomaly zone in middle crust (TLVZ)
Every method can recover the low velocity anomaly and the crust-mantle boundary (figure 3a). The misfit of MFO, SFO, and CPS330 are equal to 1.8, 3.5, and 6.3, respectively. The depth of the anomaly interfaces was accurately defined by all three methods. However, the velocity values of inverted anomaly are different from the true model. The average velocity at depth below an anomaly until the Moho of SFO and CPF330 are slightly higher than that of the true model. While MFO provides a good correlation of the inverted model with the true model. Among the three methods, SFO has the largest difference of the inverted model from true model (6.16 km/s). The SFO results show the less plausible inverted model compare to other methods. The model from SFO provides too small velocity for anomaly value. But, the mantle velocity of SFO model is much higher than that of the true values. In opposite, the predicted receiver function of MFO fit well to the data in length of ±σ (figure 4a). Both predicted receiver functions of SFO and CPS330 cannot fit the conversion phase after 5s. In all tasks, MFO provides the best quality and quantity of the results.

Broad low-velocity anomaly zone in middle crust (BLVZ)
The BLVZ model is challenging for the inversion process because the Moho depth of true model and initial model are not equal.  3b). Its predicted receiver function (figure 4b) provide the best fitting (lowest misfit) to the data. For the anomaly of the broad low-velocity zone, every inverted model can present a relative feature from the true model. However, the absolute value of the crustal velocity structure is about 0.25 km/s higher than the true model. This characteristic of the inverted model was discussed by Ammon [2]. The tradeoff between the average velocity of the model and the anomaly depth is governing the characteristic of the inverted model. The result presented from this model might suggest that neither MFO nor linearized inversion can overcome the non-uniqueness of this feature. Both methods provide a well-fit predicted receiver function, but their inverted models are insufficiently accurate.

Thin high-velocity anomaly zone in middle crust (THVZ)
Although, the anomaly of this model is opposing to the first model, TLVZ. The result has shown a similar performance ranking (figure 3c). The misfit of MFO, SFO, and CPS330 are 1.1, 1.7 and 1.2, respectively. Their predicted receiver functions (figure 4c) fit well to the synthetic data except for SFO method that waveform after 8s does not match to the synthetic data. Both MFO and CPS330 provide only a small model difference from the true model. All vital structures in the true model such as high velocity anomaly and Moho are well recovered by MFO and CPS330. In opposite, SFO shown largest model difference from the true model. Notably, the main Ps (4s) and multiple (13s) phase amplitude from the Moho of the predicted receiver function of SFO is bigger than the synthetic data. This feature is corresponding to significantly high-velocity contrast that appears at the Moho boundary of SFO inverted model. This directly confirms that the usage of multiple frequencies increases the performance of inversion. Most of the method can recover the main feature of a high velocity zone, but the broader anomaly zone of SFO effect to the flat of receiver function waveform from 8s to12s.

Broad high-velocity anomaly zone in middle crust (BHVZ)
For the BHVZ model, the performance ranking is similar to that of the third test (table 1). The misfit of MFO, SFO, and CPS330 are 1.1, 1.5 and 1.2, respectively. The high velocity anomaly is well recovered by all methods. They give an excellent performance in both model difference and data fitting. However, when the relative performances among three method were compared, the inverted model from SFO provides the highest misfit values. The recovered velocity anomaly is also broader than the true model. The predicted receiver function waveform (figure 4d) of SFO does not match to the real data after 5s.   10

The performance of Multi-frequency bands Occam's inversion
The comparison of MFO and CPS330 from four tests indicate that MFO provides a better fitting result than that from linearized inversion CPS330. Although, it must be noted that there is no user attempt to define and optimize the parameters in the inversion. The result from BLVZ model shown that MFO provides the best fitting predicted receiver function. However, its inverted model has the largest model difference from the true model. This means that the solution from Occam's inversion might fall in a local minimum of the objective function space.
Comparing the performance between SFO and MFO, all cases have shown that MFO provides a more accurate and precise inverted model in term of absolute and relative velocity structure. The TLVZ and THVZ results suggest that inversion can recover a high velocity anomaly better than low velocity anomaly. From the fourth test, BHVZ, the feature of long wavelength structure in an anomaly and sharp velocity of Moho can be recovered together by using multi-frequency data. The test also reflects the limitation of MFO. For the second method, BLVZ, the results suggest that application of multifrequency bands data still unable to solve the non-uniqueness of depth -velocity tradeoff. The missdefined Moho depth in the initial model makes all the inversion algorithm fail to construct a correct inverted model. The assumption of the initial model is still very important for the MFO algorithm. Also, this might indicate the limitation of receiver function inversion as well. The incorporate of receiver function data with other seismological data (e.g. surface wave or seismic phase travel time) still be a key to solve the problem. The integration of multi-frequencies data provides more chance to recover the sophisticated velocity model.

Conclusion
Summary, the Occam's inversion algorithm has been applied to the receiver function inverse problem. The algorithm implement in this study also allows the multiple frequency bandpass of filtered receiver function time series to be used for the inversion process. The result from synthetic tests indicate that multi-frequency bands data always provide a better result than single frequency band data. Generally, MFO provides the best fitting of predicted receiver function. However, when the initial model has a large difference from the true model, the MFO still unable to provide an accurate inverted model. The assumption of an initial model still crucial for the result of inversion. This paper has shown that the Occam's inversion can be effectively applied to the receiver function inverse problem. However, the test with real data still required to prove the performance of this algorithm.