Multi-task unscented Kalman inversion (MUKI): a derivative-free joint inversion framework and its application to joint inversion of geophysical data

Abstract

In the geophysical joint inversion, the gradient and Bayesian Markov Chain Monte Carlo (MCMC) sampling-based methods are widely used owing to their fast convergences or global optimality. However, these methods either require the computation of gradients and easily fall into local optimal solutions, or cost much time to carry out the millions of forward calculations in a huge sampling space. Different from these two methods, taking advantage of the recently developed unscented Kalman method in computational mathematics, we extend an iterative gradient-free Bayesian joint inversion framework, i.e., Multi-task unscented Kalman inversion (MUKI). In this new framework, information from various observations is incorporated, the model is iteratively updated in a derivative-free way, and a Gaussian approximation to the posterior distribution of the model parameters is obtained. We apply the MUKI to the joint inversion of receiver functions and surface wave dispersion, which is well-established and widely used to construct the crustal and upper mantle structure of the earth. Based on synthesized and real data, the tests demonstrate that MUKI can recover the model more efficiently than the gradient-based method and the Markov Chain Monte Carlo method, and it would be a promising approach to resolve the geophysical joint inversion problems.

\draftfalse\journalname

XXXXXXXXX

State Key Laboratory of Lithospheric Evolution, Institute of Geology and Geophysics, Chinese Academy of Sciences College of Earth and Planetary Sciences, University of Chinese Academy of Sciences, Beijing, China CAS Center for Excellence in Deep Earth Science, Guangzhou, China Department of Earth Science, University of Toronto, Toronto, Canada Hubei Subsurface Multi‐Scale Imaging Key Laboratory, Institute of Geophysics and Geomatics, China University of Geosciences, Wuhan, China Key Laboratory of Earth and Planetary Physics, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing, China

\correspondingauthor

Wang [email protected] \correspondingauthorChen [email protected]

{keypoints}

Multi-task unscented Kalman inversion is developed based on the unscented Kalman inversion.

Multi-task unscented Kalman inversion is a derivative-free joint inversion framework and can be applied to the joint inversion of multi-physical data.

Tests of the joint inversion of receiver functions and surface wave demonstrate MUKI’s fast convergence rate has fast-speed convergences with uncertainty estimation analysics.

Plain Language Summary

There are lots of gradient-based optimization methods and Bayesian sampling methods used in joint inversion problems. However, each of them has some disadvantages. Gradient-based methods are usually efficient in fast convergence, but they easily fall into local minima; Bayesian sampling methods can help to seek global solutions but are computationally intensive. To take their advantages simultaneously, here we propose a new joint inversion framework, Multi-task unscented Kalman inversion (MUKI), based on the recently developed unscented Kalman inversion method. In this new framework, a gradient-free, iterative manner is adopted to yield a fast convergence, and a Bayesian approach is used for uncertainty analysis. Joint inversion of receiver functions and surface wave dispersion is well-established in geophysics and has been successfully applied with gradient-based methods or Bayesian sampling methods to reconstruct the structure of the earth’s crust and upper mantle. Thus, we use the joint inversion of receiver functions and surface wave as a benchmark for our newly developed MUKI. The tests of synthesized and real data prove that MUKI is more robust, accurate, and efficient compared to the gradient and Bayesian sampling-based methods.

1 Introduction

Geophysical inversion is an optimization problem that uses mathematical methods to estimate geological model parameters from a series of observed geophysical data [Fichtner (\APACyear2021), Stuart (\APACyear2010), Tarantola (\APACyear2005)]. Independent inversions of a single property are ineffective when the resolution of a single data set is limited. Inverting one geophysical data set with other supplemental geophysical data sets simultaneously, i.e., joint inversion, can improve inversion results by utilizing different sensitivities of different data sets to specific parameters and can mitigate the nonuniqueness of inversion problems compared with only using a single data set.

Joint inversion can be divided into two groups, either multiple data sets that are sensitive to the same physical properties or different data sets that are responses to the different physical properties of the same geological structure [Moorkamp \BOthers. (\APACyear2007)]. These two groups can both be solved by gradient-based methods, such as the damped least squares method [Gallardo (\APACyear2007), Julià \BOthers. (\APACyear2000)]. Although gradient-based methods can be effective, these methods easily suffer from local minimum except for a good initial model. To obtain geological realistic models, these methods usually adopt some damping and smooth regularization, which may reduce accuracy or resolution. The recently developed Bayesian-based MCMC method can not only obtain a global optimal solution but also obtain uncertainties of results. However, the MCMC method is costly because millions of forward calculations are required to thoroughly sample the posterior probability density function.

Overall, combining different model parameters with observed data sets for an effective joint inversion method remains a challenge. To solve this problem, we resort to a new inversion framework, i.e., Kalman filter-based (KF) inversion. KF was first proposed to describe a recursive solution to the discrete-data linear filtering problem in 1960 [Kalman (\APACyear1960)]. Since then, KF has been a subject of extensive research and application, particularly in the area of signal processing [Dean (\APACyear1986), Jwo \BBA Lai (\APACyear2008), Plett (\APACyear2006)].

Two Kalman-based inversion methods can be applied to nonlinear systems, ensemble Kalman inversion (EnKI) and unscented Kalman inversion (UKI), which can be seen as applications of the corresponding filters in the optimization domain. KF itself is intended for linear systems. When the system is nonlinear, this method can be extended to the extended Kalman filter (ExKF) by linearizing the nonlinear system. However, the ExKF method is only reliable for weakly nonlinear problems. To resolve major problems of the ExKF for parameter estimation for strongly nonlinear models and large state spaces, \citeAEvensen1994 proposed a Monte Carlo-based KF, i.e., ensemble Kalman filter. By using a collection of hundreds to thousands of state vectors to evaluate the system, this novel filter works well in both linear and nonlinear systems [Evensen (\APACyear2009), Iglesias \BOthers. (\APACyear2012)], and there are several optimization implementations in recent years [Aleardi \BOthers. (\APACyear2021), Conjard \BBA Grana (\APACyear2021), Muir \BBA Tsai (\APACyear2019), Wang \BOthers. (\APACyear2019)].

Analogous to EnKI, UKI is also gradient-free but adopts some integration points (which are the deterministic sample points to capture the mean and covariance of the Gaussian distribution) to evaluate a nonlinear system. By using the so-called unscented transformation [Julier \BOthers. (\APACyear1995)], UKI has been shown to produce superior results to EnKI on a variety of nonlinear inverse problems [Huang \BOthers. (\APACyear2022)].

Here we extend UKI to multi-task UKI (MUKI) by introducing a generalized covariance function across different data sets and model parameters. As a benchmark, we apply the MUKI to the joint inversion of receiver functions and surface wave dispersion, which is a powerful tool to invert the crustal and upper mantle structure of the earth, and the gradient-based method and MCMC method have been successfully intro[Julià \BOthers. (\APACyear2003), Sambridge \BOthers. (\APACyear2013)]. The remainder of this paper will introduce and evaluate our developed MUKI method comprehensively against these classic methods.

2 Method

2.1 Gaussian Approximation Algorithm

In general, geophysical inverse problems can be formulated as follows

d=G(m)+\eta,

(1)

where $G:\mathbb{R}^{N_{m}}\rightarrow\mathbb{R}^{N_{d}}$ is a nonlinear function (especially, a forward map in RF or SWD) that denotes a forward operator. In inversion, the data $d$ is available, thus the inverse problem is to recover parameter vector $m\in\mathbb{R}^{N_{m}}$ from $d\in\mathbb{R}^{N_{y}}$ , with the observed noise $\eta$ drawn from a Gaussian distribution $\mathcal{N}(0,\Sigma_{\eta})$ .

In the Bayesian viewpoint, the inverse process aims to maximize a posterior estimation for $m$ , which can be written as

p(m|D)\propto\exp(-\Phi(m))p_{0}(m)\hskip 7.22743pt\text{ with }\hskip 7.22743pt\Phi(m)=\frac{1}{2}\left\|\Sigma_{\eta}^{-\frac{1}{2}}(d-G(m))\right\|^{2}.

(2)

Generally, the exact posterior distribution $p(m|D)$ , is intractable to compute. In this paper, we view the process of maximizing a posterior estimation as a stochastic dynamic for the parameter, so that we can employ KF to estimate the parameter for a given observation [Huang \BBA Huang (\APACyear2021)].

Consider the stochastic dynamic system:

	$\displaystyle evolution:m_{n+1}=m_{n}+\omega_{n+1},\omega_{n+1}\sim\mathcal{N}(0,\Sigma_{\omega}),$		(3)
	$\displaystyle observation:d_{n+1}=G(m_{n+1})+\nu_{n+1},\nu_{n+1}\sim\mathcal{N}(0,\Sigma_{\nu}),$		(4)

where the artificial evolution error covariance $\Sigma_{\omega}$ and the artificial observation error covariance $\Sigma_{\nu}$ are both positive definite matrices. We use the Gaussian Approximation Algorithm (GAA) to obtain $m_{n}$ from $D_{n}={d_{1},d_{2},...d_{n}}$ , where $d_{n}$ is the observation $d$ at $n$ iteration. Note that when we apply KF in optimization, the $d_{n}$ is identical to $d_{obs}$ [Huang \BOthers. (\APACyear2022)].

The GAA, originally in geostatistics, also known as Gaussian process regression [Chiles \BBA Delfiner (\APACyear2012)], is a method of interpolation that maps Gaussian distribution to Gaussian distribution during a Gaussian process and leads to an insight into the Kalman methodology [Huang \BOthers. (\APACyear2022)]. This algorithm starts with a Gaussian approximation to the posterior distribution $p_{0}\sim\mathcal{N}(\mu_{0},C_{0})$ , $\mu_{0}$ and $C_{0}$ are the mean and covariance of the Gaussian distribution $p_{0}$ , respectively. Then successive estimates, $p_{1},p_{2},...p_{n}$ can be iteratively computed by repeatedly applying the following two steps, $p_{n}\mapsto\hat{p}_{n+1}$ and then $\hat{p}_{n+1}\mapsto p_{n+1}$ . In the first step (analogue to prediction step in Kalman Filter), the predicted model parameters vector $\hat{p}_{n+1}\sim\mathcal{N}\left(\tilde{\mu}_{n+1},\tilde{C}_{n+1}\right)$ is also Gaussian and satisfies

	$\displaystyle\tilde{\mu}_{n+1}$	$\displaystyle=$	$\displaystyle\mu_{n},$		(5)
	$\displaystyle\tilde{C}_{n+1}$	$\displaystyle=$	$\displaystyle C_{n}+\Sigma_{\omega},$		(6)

where $\Sigma_{\omega}$ are observation error covariances, and $\tilde{.}$ denotes the updated parameters. Then in the second step (correction step), we use the mean and covariance matrix to represent the joint distribution of $\left\{m_{n+1},d_{n+1}\mid D_{n}\right\}$ , which is a multivariate Gaussian distribution

\displaystyle\mathcal{N}\left(\left[\begin{array}[]{lcl}\tilde{\mu}_{n+1}\\ \tilde{d}_{n+1}\end{array}\right],\left[\begin{array}[]{lcl}\tilde{C}_{n+1}&\tilde{C}^{md}_{n+1}\\ \tilde{C}^{dm}_{n+1}&\tilde{C}^{dd}_{n+1}\end{array}\right]\right),

(11)

where

$\displaystyle\tilde{d}_{n+1}$	$\displaystyle=$	$\displaystyle\mathbb{E}[d_{n+1}\|D_{n}]=\mathbb{E}[G(m_{n+1})\|D_{n}],$	(12)
$\displaystyle\tilde{C}^{md}_{n+1}$	$\displaystyle=$	$\displaystyle Cov[m_{n+1},d_{n+1}\|D_{n}]=Cov[m_{n+1},G(m_{n+1})\|D_{n}],$	(13)
$\displaystyle\tilde{C}^{dd}_{n+1}$	$\displaystyle=$	$\displaystyle Cov[d_{n+1}\|D_{n}]=Cov[G(m_{n+1})\|D_{n}]+\Sigma_{\nu}.$	(14)

$\Sigma_{\nu}$ is the artificial observation error covariance.

Conditioning the joint Gaussian distribution in equation (11) to obtain p( $m_{n+1}|D_{n+1}=d_{n+1})$ [Bishop (\APACyear2006)], specific operations can also be found in Supplementary Text S1.

	$\displaystyle\mu_{n+1}$	$\displaystyle=$	$\displaystyle\tilde{\mu}_{n+1}+\tilde{C}^{md}_{n+1}(\tilde{C}^{dd}_{n+1})^{-1}(d_{n+1}-\tilde{d}_{n+1}),$		(15)
	$\displaystyle C_{n+1}$	$\displaystyle=$	$\displaystyle\tilde{C}_{n+1}-\tilde{C}^{md}_{n+1}(\tilde{C}^{dd}_{n+1})^{-1}\tilde{C}^{md}_{n+1}{}^{T}.$		(16)

We can see that GAA not only provides the expected value $\mu_{n+1}$ for the given observation $d_{n+1}$ , but also offers uncertainty through the variance $C_{n+1}$ . From an optimization point of view, the observation $D_{n}$ is not changed, it is identical to $d_{obs}$ at each iteration. By implementing approximations to equations (12) to (14), one iteration of the Kalman inversion will be established.

In this research, the parameters in equations (6), (14) are chosen as $\Sigma_{\omega}=C_{n}$ and $\Sigma_{v}=2\Sigma_{\eta}$ , which guarantees that the algorithm can obtain an accurate Gaussian estimation to the posterior probability with the converged mean and covariance [Huang \BBA Huang (\APACyear2021)].

2.2 Unscented Kalman Inversion (UKI)

In the Kalman method, a vital operation is how to measure the statistical properties of the system states, i.e. approximating the equations (13),(14). The UKI method approximates the integrals by deterministic quadrature rules, which are generally called unscented transform [Julier \BOthers. (\APACyear1995)]. The approach is also described in Supplementary Text S2. To evaluate the integration, 2 $N_{m}+1$ ( $N_{m}$ is the dimension of the vector $m$ ) integration points have been selected deterministically according to the posterior distribution $\mathcal{N}(\mu_{n+1},C_{n+1})$

	$\displaystyle m_{n+1}^{j}$	$\displaystyle=$	$\displaystyle\mu_{n+1}+c_{j}[\sqrt{C_{n+1}}]_{j}\hskip 14.45377pt1\leq j\leq N_{m},$		(17)
	$\displaystyle m_{n+1}^{j+N_{m}}$	$\displaystyle=$	$\displaystyle\mu_{n+1}-c_{j}[\sqrt{C_{n+1}}]_{j}\hskip 14.45377pt1\leq j\leq N_{m}.$		(18)

where $[\sqrt{C_{n+1}}]_{j}$ is the j-th column of the Cholesky factor of $C_{n+1}$ in the $n+1$ step, according to the parameter setting in unscented transform, the coefficient $c_{j}=a\sqrt{N_{m}}$ and the hyperparameters $a=min\{\sqrt{\frac{4}{N_{m}}},1\}$ . The central integration point $m_{n+1}^{0}=\mu_{n+1}$ .

The equations (12) to (14) can be defined as

$\displaystyle\tilde{d}_{n+1}$	$\displaystyle=$	$\displaystyle G(m_{n+1}^{0}),$	(19)
$\displaystyle\tilde{C}^{md}_{n+1}$	$\displaystyle=$	$\displaystyle\sum_{j=1}^{2N_{m}}w_{j}^{c}(m_{n+1}^{j}-\mu_{n+1})(m_{n+1}^{j}-\mathbb{E}G(\mu_{n+1}))^{T},$	(20)
$\displaystyle\tilde{C}^{dd}_{n+1}$	$\displaystyle=$	$\displaystyle\sum_{j=1}^{2N_{m}}w_{j}^{c}(G(m_{n+1}^{j})-\mathbb{E}G(\mu_{n+1}))(G(m_{n+1}^{j})-\mathbb{E}G(\mu_{n+1}))^{T},$	(21)

where quadrature weights $w_{j}^{c}=\frac{1}{2a^{2}N_{m}}$ . Details on this approach can be found in \citeAHuang2021a and the supplementary Text S2.

Refer to caption — Figure 1: Schematical map of the MUKI framework. This diagram only shows the process of one model parameter update, the loop termination condition is not drawn and can be set according to the error function.

2.3 Multi-task unscented Kalman inversion

Inspired by the multi-task learning approach in machine learning problems [Bonilla \BOthers. (\APACyear2007)], in this study, we combine the observation data $d_{A},d_{B},d_{C},...$ (i.e., the multi-responses of models), into a generalized $d$ , so that different observations (we only use $d_{A},d_{B},d_{C}$ to demonstratively represent multi-responses of models in the remaining part of the paper) can share covariance. The generalized $d$ is set as

\displaystyle d=[\alpha d_{A1},...,\alpha d_{A_{m}},\beta d_{B_{1}},...,\beta d_{B_{n}},\gamma d_{C1},...,\gamma d_{Co},...]^{T},

(22)

where $\alpha,\beta,\gamma,...$ are custom weight factors according to the relative quality of data and weight.

We can summarize the process of the MUKI inversion in Figure 1. The process shown in Figure 1 can be divided into two steps, the prediction step (grey background) and the correction step (white background), which correspond to the equations (5) to (6) and (13)-(16), respectively. In the correction step, there are five subsections, including the predicted model parameters, the forward data, the observed data, the two covariance matrices, and updated equations, i.e., equations (15) and (16). These white circles represent random variables, i.e. four model parameters and four observation parameters. It is strait forward to multiple parameters. The predicted model parameters will be predicted by the prediction step, the synthetic data is the theoretical data that will be generated by the forward computation. Note that $d_{A}$ can be generated by $m_{\alpha}$ , or by $m_{\alpha}$ , $m_{\beta}$ simultaneously, which relies on whether $d_{A}$ depends on $m_{\beta}$ . Our model mainly shares correlations through the covariance matrix, which will be generated by the unscented transformation. If $C^{m_{\alpha}d_{B}}=0$ , it means that the parameters $m_{a}lpha$ and $d_{B}$ are not correlated and if $C^{m_{\alpha}d_{B}}$ , $C^{m_{\beta}d_{A}}$ , $C^{d_{AB}}$ are all 0, the joint inversion will degenerate to an individual inversion of the two physical processes.

Although the Kalman inversion does not need to optimize (minimize) an objective function, we use the squared error cost function to purely measure the convergence of the predicted model vector $m$

\displaystyle\Phi(m)=\frac{1}{2}\left\|\Sigma_{\eta}^{-\frac{1}{2}}(d-G(\mu))\right\|^{2},

(23)

where observed noise $\eta=[\alpha\eta_{A},\beta\eta_{B},\gamma\eta_{C},...]^{T}$ . Given a custom weight is imposed on the data, thus the objective function is defined by the linear combination of weighted misfits, thus taking the form:

\displaystyle\Phi(m)=\alpha\Phi_{A}+\beta\Phi_{B}+\gamma\phi_{C}+....

(24)

3 Results

Joint inversion of receiver functions (RFs) and surface wave dispersion (SWD) is well-established in geophysics and has been successfully applied with gradient-based methods and Bayesian sampling methods to recover the structure of the earth’s crust and upper mantle. Thus, we use the joint inversion of RFs (only use radial receiver function in this paper) and SWD as a benchmark for our newly developed MUKI to verify its effectiveness. Although it is a general inversion framework for joint inversion, we set the number of observations as 2, i.e., $d_{A}$ and $d_{B}$ are RF and SWD data respectively. The model is set up as a horizontally layered model beneath the receivers with the lowest layer a half space. The model consists of 4 parameters, i.e., $m=[v_{p},v_{s},\rho,h]$ , where $v_{p},v_{s},\rho$ and $h$ are respectively the P wave velocity, S wave velocity, density and thickness at each given layer. In the joint inversion, $v_{p}$ and $\rho$ can be determined by $v_{s}$ with empirical relations [Brocher (\APACyear2005)], which can also be found in Supplementary Text S3. Due to the weak sensitivity of P wave velocity and density to data, we just invert the S wave velocity and thickness of the model.

We set thickness and S wave velocity value in each layer are both variables so that their posterior distributions can be estimated through equations 15 and 16. In the joint inversion of RFs and SWD, the mean of the model parameter is initialized as $m=[\mu_{v}s,\mu_{t}hickness]^{T}$ .

We demonstrate the effectiveness and efficiency of the MUKI using several tests to invert synthetic and real data of RF(s) and SWD. Throughout all applications, we focus on MUKI. Some comparisons with the UKI, MCMC [Bodin, Sambridge, Rawlinson\BCBL \BBA Arroucau (\APACyear2012)], and gradient-based methods for the joint inversion of RFs and SWD (CPS,\citeAHerrmann2013) are also considered.

3.1 Synthetic Example

We use two numerical velocity models to evaluate our algorithm, including i) Layered model, which consists of 8 horizontal layers with a low S wave velocity layer in the crust and a sharp velocity gradient at the Moho. The model is modified from the 6 layers model referring to \citeABodin2012a, ii) GAr1 model, which can be viewed as a smoother version of model i [Sambridge (\APACyear1999)]. We use the matrix-propagation method to synthesize the seismograms and then obtain the RFs and SWD data [Haskell (\APACyear1953), Herrmann (\APACyear2013)]. Specially, we use the water-level deconvolution method to calculate the RF with a water-level factor of 0.001. Following the work of \citeABodin2012, we assume the RF or SWD noise is drawn from the multivariate Gaussian distribution and can be parameterized with two parameters. For SWD, the noise is generated by a diagonal covariance matrix (i.e., n $\sigma_{true}^{SWD}=0.012$ and $r_{true}^{SWD}=0$ ); while RF, the noise is generated with an exponential correlation law with values $\sigma_{true}^{RF}=0.005$ and $r_{true}^{RF}=0.92$ .

To initialize MTUKI, we set that the initial model consists of 25 layers, and the means of the thickness of the first seventh layers are all 2km, so the $\mu_{thk}$ can be set as $\mu_{thk}=[2,2,2,2,2,2,2,3,...,3]$ . Considering the covariance of the parameter can be iterated updated in the inversion process as equation (16), we choose $\mu_{0}=[\mu_{vs},\mu_{thk}],C_{0}=0.001\mathbb{I}$ and initial MUKI at $m_{0}\sim\mathcal{N}[\mu_{0},C_{0}]$ , where $\mu_{0}$ is the prior mean and $C_{0}$ is an uninformative prior covariance matrix. Note that $\mu_{0}=[\mu_{vs},\mu_{thk}]$ , and the setting of $\mu_{vs}$ will be discussed later.

Figure 1 and Figure 2 show the inversion processes in the 8-layer model. In Figure 1, it can be seen that the algorithm can adaptively adjust the model parameters (i.e.velocity and thickness of each layer) (Figure 1a) and fits the observations well (Figure 1b,c). The convergence of the parameter vector is shown in Figure 2d, Figure 2 shows the profiles of the mean and standard deviation (std) of Vs, it can be observed that this method can converge to an ”optimal” solution with only $O(10)$ iterations (while CPS need 25 iterations). Compared with the MCMC method (Figure 2b) and gradient method (Figure 2c), our method can obtain an accurate model at least above 60 km with a shorter CPU time (6.01 s/6866.91 s/20 s) on the same computing environment (Intel Core i7-10875H). One detail that must be considered is that our method considers both the thickness and S wave velocity are variables, so the depth of our inversion result is usually different from the initial model setting. Due to the CPS method cannot access the error information directly, we follow the work of \citeAliGeodynamicProcessesContinental2020 to regenerate theoretical RF and SWD by adding Gaussian random noise as described above, then we obtain the 110 data sets (one RF and one SWD) and then performed the joint inversion respectively to obtain 110 velocity models. After the statistical analysis, we show the results in 2c. Since CPS performs the joint inversion using a damped least squares algorithm [Julià \BOthers. (\APACyear2000)] that contains a priori smoothness constraint on the velocity of adjacent layers, thus the results may be smoother than the reference model.

The comparisons between the joint inversion result with the results of individual RF and SWD inversions are shown in Figure S1, which shows that the joint inversion can obtain the best parameter estimation against the individual inversions of SWD and RF.

Then we design several tests to investigate the influence of the initial Vs on our result (Figure S2). In these tests, the initial $\mu_{thk}$ are identical but $\mu_{vs}$ are different. This numerical experiment demonstrates that discontinuities with sharp speed changes (i.e., Moho interface) can be effectively recovered with different initial $\mu_{vs}$ .

3.2 Field Case

Our algorithm is further tested using a real data set from the KIGAM station collected by Seoul National University, which is located at $(35.00^{\circ}N,126.249^{\circ}E)$ . The dispersion measurements and receiver functions at station KIGAM have been obtained and performed well [Herrmann (\APACyear2013)]. We performed experiments to compare the MUKI method with the other two methods. The RF was calculated using the time domain iterative deconvolution procedure [Ligorría \BBA Ammon (\APACyear1999)], which is more stable in the presence of noisy data than the frequency domain method. The RF has a Gaussian factor alpha=2.5 and the period of SWD (Phase velocity) varies from 10 s to 37 s. The initial mean of thickness and S wave velocity are shown in Figure 4a, which is the same as the setting in CPS.

The joint inversion results for the station KIGAM are shown in Figure 4. Figure 4a shows that it converges after twenty iterations, which are drawn from light blue line to dark blue line, and can invert the distribution of thickness and S wave velocity in each iteration. Since we only use the surface wave dispersion of less than 40 s, the uncertainty of the inversion results will increase rapidly with depth when the depth exceeds 80 km. To illustrate this, we have plotted in Figure 4b the uncertainty for the S wave velocity of each layer. As the depth increases, the maximum probability density of S wave velocity at each depth will decrease, which means an increase in uncertainty. As RF has strong non-uniqueness for the shallow S wave velocity, leading to large uncertainties. The solutions obtained by the MUKI (red dashed line), CPS (green line), and MCMC (blue line) are shown in Figure 4c, respectively. The probability density of the posterior distribution of Vs is plotted as the base map. In Figures 4d and 4e, it can be found that the structure recovered by MUKI can fit the observed RF and SWD best, which can also be confirmed by the optimization error in Figures 4f, and 4g. These results verify that the model recovered by our method in Figure 4b is more reliable, and the discontinuous interface around 32-34 km (Moho layer [Chang \BBA Baag (\APACyear2007)]) can be recovered well.

4 Conclusions

In this study, we develop a new derivative-free joint inversion framework and evaluate it by joint inversion of RF and SWD. In both synthetic and real data set tests, the newly developed inversion framework demonstrates its powerful ability in terms of at least three aspects: a) Due to it is derivative-free inversion framework, it honors good flexibility in inverting multi-physical geophysical data; b) It can effectively obtain uncertainty of the solution by Gaussian approximation; c) It has a fast convergence rate compared to traditional methods.

5 Open Research

RF data and the SWD data from the station KIGAM can be accessed through the CPS tutorial (http://www.eas.slu.edu/eqc/eqc_cps/TUTORIAL/STRUCT/index.html).

Acknowledgements.

This research was jointly supported by the National Key R

\&

D Program of China (grant 2016YFC0600402), the Strategic Priority Research Program (B) of the Chinese Academy of Sciences (grant XDB18000000), and the National Natural Science Foundation of China (grants 41374063 and 41874065). We appreciate Daniel Z. Huang for his outstanding work about UKI. Thanks to Ran You for his great helps with UKF. We are grateful to R.B. Herrmann and Jennifer Dreiling for providing the CPS software and MCMC code used in this study. Thanks to Sicheng Zuo for revising the manuscript. Thanks to Wentao Li and Yifan Lu for the valuable discussion about the manuscript. Thanks to colleagues in our lab for helpful discussions on the manuscript.

References

Aleardi \BOthers. (\APACyear2021) \APACinsertmetastarAleardi2021{APACrefauthors}Aleardi, M., Vinciguerra, A.\BCBL \BBA Hojat, A. \APACrefYearMonthDay2021. \BBOQ\APACrefatitleEnsemble-Based Electrical Resistivity Tomography with Data and Model Space Compression Ensemble-based electrical resistivity tomography with data and model space compression.\BBCQ \APACjournalVolNumPagesPure and Applied Geophysics17851781–1803. {APACrefDOI} 10.1007/s00024-021-02730-1 \PrintBackRefs\CurrentBib
Bishop (\APACyear2006) \APACinsertmetastarBishop2006{APACrefauthors}Bishop, C\BPBIM. \APACrefYear2006. \APACrefbtitlePattern Recognition and Machine Learning Pattern recognition and machine learning. \APACaddressPublisherSpringer. \PrintBackRefs\CurrentBib
Bodin, Sambridge, Rawlinson\BCBL \BBA Arroucau (\APACyear2012) \APACinsertmetastarBodin2012{APACrefauthors}Bodin, T., Sambridge, M., Rawlinson, N.\BCBL \BBA Arroucau, P. \APACrefYearMonthDay2012. \BBOQ\APACrefatitleTransdimensional Tomography with Unknown Data Noise Transdimensional tomography with unknown data noise.\BBCQ \APACjournalVolNumPagesGeophysical Journal International18931536–1556. {APACrefDOI} 10.1111/j.1365-246X.2012.05414.x \PrintBackRefs\CurrentBib
Bodin, Sambridge, Tkalčić\BCBL \BOthers. (\APACyear2012) \APACinsertmetastarBodin2012a{APACrefauthors}Bodin, T., Sambridge, M., Tkalčić, H., Arroucau, P., Gallagher, K.\BCBL \BBA Rawlinson, N. \APACrefYearMonthDay2012. \BBOQ\APACrefatitleTransdimensional Inversion of Receiver Functions and Surface Wave Dispersion: TRANSDIMENSIONAL INVERSION OF RF AND SWD Transdimensional inversion of receiver functions and surface wave dispersion: Transdimensional inversion of rf and swd.\BBCQ \APACjournalVolNumPagesJournal of Geophysical Research: Solid Earth117B2B02301. {APACrefDOI} 10.1029/2011JB008560 \PrintBackRefs\CurrentBib
Bonilla \BOthers. (\APACyear2007) \APACinsertmetastarBonilla2007{APACrefauthors}Bonilla, E\BPBIV., Chai, K.\BCBL \BBA Williams, C. \APACrefYearMonthDay2007. \BBOQ\APACrefatitleMulti-Task Gaussian Process Prediction Multi-task gaussian process prediction.\BBCQ \BIn \APACrefbtitleAdvances in Neural Information Processing Systems Advances in neural information processing systems (\BVOL 20). \PrintBackRefs\CurrentBib
Brocher (\APACyear2005) \APACinsertmetastarBrocher2005{APACrefauthors}Brocher, T\BPBIM. \APACrefYearMonthDay2005. \BBOQ\APACrefatitleEmpirical Relations between Elastic Wavespeeds and Density in the Earth’s Crust Empirical relations between elastic wavespeeds and density in the earth’s crust.\BBCQ \APACjournalVolNumPagesBulletin of the Seismological Society of America9562081–2092. {APACrefDOI} 10/cn3qwb \PrintBackRefs\CurrentBib
Chang \BBA Baag (\APACyear2007) \APACinsertmetastarChang2007{APACrefauthors}Chang, S\BHBIJ.\BCBT \BBA Baag, C\BHBIE. \APACrefYearMonthDay2007. \BBOQ\APACrefatitleMoho Depth and Crustal VP/VS Variation in Southern Korea from Teleseismic Receiver Functions: Implication for Tectonic Affinity between the Korean Peninsula and China Moho depth and crustal vp/vs variation in southern korea from teleseismic receiver functions: Implication for tectonic affinity between the korean peninsula and china.\BBCQ \APACjournalVolNumPagesBulletin of the Seismological Society of America9751621–1631. {APACrefDOI} 10.1785/0120050264 \PrintBackRefs\CurrentBib
Chiles \BBA Delfiner (\APACyear2012) \APACinsertmetastarChiles2012{APACrefauthors}Chiles, J\BHBIP.\BCBT \BBA Delfiner, P. \APACrefYear2012. \APACrefbtitleGeostatistics: Modeling Spatial Uncertainty Geostatistics: Modeling spatial uncertainty (\PrintOrdinal2nd ed \BEd). \APACaddressPublisherWiley. \PrintBackRefs\CurrentBib
Conjard \BBA Grana (\APACyear2021) \APACinsertmetastarConjard2021{APACrefauthors}Conjard, M.\BCBT \BBA Grana, D. \APACrefYearMonthDay2021. \BBOQ\APACrefatitleEnsemble-Based Seismic and Production Data Assimilation Using Selection Kalman Model Ensemble-Based Seismic and Production Data Assimilation Using Selection Kalman Model.\BBCQ \APACjournalVolNumPagesMathematical Geosciences5371445–1468. {APACrefDOI} 10.1007/s11004-021-09940-2 \PrintBackRefs\CurrentBib
Dean (\APACyear1986) \APACinsertmetastarDean1986{APACrefauthors}Dean, G\BPBIC. \APACrefYearMonthDay1986. \BBOQ\APACrefatitleAn Introduction to Kalman Filters An introduction to kalman filters.\BBCQ \APACjournalVolNumPagesMeasurement and Control1969–73. {APACrefDOI} 10.1177/002029408601900204 \PrintBackRefs\CurrentBib
Evensen (\APACyear1994) \APACinsertmetastarEvensen1994{APACrefauthors}Evensen, G. \APACrefYearMonthDay1994. \BBOQ\APACrefatitleSequential Data Assimilation with a Nonlinear Quasi-Geostrophic Model Using Monte Carlo Methods to Forecast Error Statistics Sequential data assimilation with a nonlinear quasi-geostrophic model using monte carlo methods to forecast error statistics.\BBCQ \APACjournalVolNumPagesJournal of Geophysical Research: Oceans99C510143–10162. {APACrefDOI} 10.1029/94JC00572 \PrintBackRefs\CurrentBib
Evensen (\APACyear2009) \APACinsertmetastarEvensen2009{APACrefauthors}Evensen, G. \APACrefYear2009. \APACrefbtitleData Assimilation Data assimilation. \APACaddressPublisherSpringer Berlin Heidelberg. {APACrefDOI} 10.1007/978-3-642-03711-5 \PrintBackRefs\CurrentBib
Fichtner (\APACyear2021) \APACinsertmetastarFichtner2021{APACrefauthors}Fichtner, A. \APACrefYear2021. \APACrefbtitleLecture Notes on Inverse Theory Lecture notes on inverse theory [preprint]. {APACrefDOI} 10.33774/coe-2021-qpq2j \PrintBackRefs\CurrentBib
Gallardo (\APACyear2007) \APACinsertmetastarGallardo2007{APACrefauthors}Gallardo, L\BPBIA. \APACrefYearMonthDay2007. \BBOQ\APACrefatitleMultiple Cross-Gradient Joint Inversion for Geospectral Imaging Multiple cross-gradient joint inversion for geospectral imaging.\BBCQ \APACjournalVolNumPagesGeophysical Research Letters3419. {APACrefDOI} 10.1029/2007GL030409 \PrintBackRefs\CurrentBib
Haskell (\APACyear1953) \APACinsertmetastarHaskell1953{APACrefauthors}Haskell, N\BPBIA. \APACrefYearMonthDay1953. \BBOQ\APACrefatitleThe Dispersion of Surface Waves on Multilayered Media* The dispersion of surface waves on multilayered media*.\BBCQ \APACjournalVolNumPagesBulletin of the Seismological Society of America43117–34. \PrintBackRefs\CurrentBib
Herrmann (\APACyear2013) \APACinsertmetastarHerrmann2013{APACrefauthors}Herrmann, R\BPBIB. \APACrefYearMonthDay2013. \BBOQ\APACrefatitleComputer Programs in Seismology: An Evolving Tool for Instruction and Research Computer programs in seismology: An evolving tool for instruction and research.\BBCQ \APACjournalVolNumPagesSeismological Research Letters8461081–1088. \PrintBackRefs\CurrentBib
Huang \BBA Huang (\APACyear2021) \APACinsertmetastarHuang2021a{APACrefauthors}Huang, D\BPBIZ.\BCBT \BBA Huang, J. \APACrefYearMonthDay2021. \APACrefbtitleUnscented Kalman Inversion: Efficient Gaussian Approximation to the Posterior Distribution. Unscented kalman inversion: Efficient gaussian approximation to the posterior distribution. \PrintBackRefs\CurrentBib
Huang \BOthers. (\APACyear2022) \APACinsertmetastarHuang2022{APACrefauthors}Huang, D\BPBIZ., Schneider, T.\BCBL \BBA Stuart, A\BPBIM. \APACrefYearMonthDay2022. \BBOQ\APACrefatitleIterated Kalman Methodology for Inverse Problems Iterated kalman methodology for inverse problems.\BBCQ \APACjournalVolNumPagesJournal of Computational Physics463111262. {APACrefDOI} 10.1016/j.jcp.2022.111262 \PrintBackRefs\CurrentBib
Iglesias \BOthers. (\APACyear2012) \APACinsertmetastarIglesias2012{APACrefauthors}Iglesias, M\BPBIA., Law, K\BPBIJ\BPBIH.\BCBL \BBA Stuart, A\BPBIM. \APACrefYearMonthDay2012. \BBOQ\APACrefatitleEnsemble Kalman Methods for Inverse Problems Ensemble kalman methods for inverse problems.\BBCQ \APACjournalVolNumPagesInverse Problems29045001. {APACrefDOI} 10.1088/0266-5611/29/4/045001 \PrintBackRefs\CurrentBib
Julier \BOthers. (\APACyear1995) \APACinsertmetastarJulier1995{APACrefauthors}Julier, S., Uhlmann, J.\BCBL \BBA Durrant-Whyte, H. \APACrefYearMonthDay1995. \BBOQ\APACrefatitleA New Approach for Filtering Nonlinear Systems A new approach for filtering nonlinear systems.\BBCQ \BIn \APACrefbtitleProceedings of 1995 American Control Conference - ACC’95 Proceedings of 1995 american control conference - acc’95 (\BVOL 3, \BPG 1628-1632 vol.3). {APACrefDOI} 10.1109/ACC.1995.529783 \PrintBackRefs\CurrentBib
Julià \BOthers. (\APACyear2003) \APACinsertmetastarJulia2003{APACrefauthors}Julià, J., Ammon, C\BPBIJ.\BCBL \BBA Herrmann, R\BPBIB. \APACrefYearMonthDay2003. \BBOQ\APACrefatitleLithospheric Structure of the Arabian Shield from the Joint Inversion of Receiver Functions and Surface-Wave Group Velocities Lithospheric structure of the Arabian Shield from the joint inversion of receiver functions and surface-wave group velocities.\BBCQ \APACjournalVolNumPagesTectonophysics37111–21. {APACrefDOI} 10.1016/S0040-1951(03)00196-3 \PrintBackRefs\CurrentBib
Julià \BOthers. (\APACyear2000) \APACinsertmetastarJulia2000{APACrefauthors}Julià, J., Ammon, C\BPBIJ., Herrmann, R\BPBIB.\BCBL \BBA Correig, A\BPBIM. \APACrefYearMonthDay2000. \BBOQ\APACrefatitleJoint Inversion of Receiver Function and Surface Wave Dispersion Observations Joint inversion of receiver function and surface wave dispersion observations.\BBCQ \APACjournalVolNumPagesGeophysical Journal International143199–112. {APACrefDOI} 10.1046/j.1365-246x.2000.00217.x \PrintBackRefs\CurrentBib
Jwo \BBA Lai (\APACyear2008) \APACinsertmetastarJwo2008{APACrefauthors}Jwo, D\BHBIJ.\BCBT \BBA Lai, C\BHBIN. \APACrefYearMonthDay2008. \BBOQ\APACrefatitleUnscented Kalman Filter with Nonlinear Dynamic Process Modeling for GPS Navigation Unscented kalman filter with nonlinear dynamic process modeling for gps navigation.\BBCQ \APACjournalVolNumPagesGPS Solutions124249–260. {APACrefDOI} 10.1007/s10291-007-0081-9 \PrintBackRefs\CurrentBib
Kalman (\APACyear1960) \APACinsertmetastarKalman1960{APACrefauthors}Kalman, R\BPBIE. \APACrefYearMonthDay1960. \BBOQ\APACrefatitleA New Approach to Linear Filtering and Prediction Problems A new approach to linear filtering and prediction problems.\BBCQ \APACjournalVolNumPagesJournal of Basic Engineering82135–45. {APACrefDOI} 10.1115/1.3662552 \PrintBackRefs\CurrentBib
Li \BOthers. (\APACyear2020) \APACinsertmetastarliGeodynamicProcessesContinental2020{APACrefauthors}Li, W., Chen, Y., Tan, P.\BCBL \BBA Yuan, X. \APACrefYearMonthDay2020. \BBOQ\APACrefatitleGeodynamic processes of the continental deep subduction: Constraints from the fine crustal structure beneath the Pamir plateau Geodynamic processes of the continental deep subduction: Constraints from the fine crustal structure beneath the pamir plateau.\BBCQ \APACjournalVolNumPagesScience China Earth Science63649–661. {APACrefDOI} 10.1007/s11430-019-9587-3 \PrintBackRefs\CurrentBib
Ligorría \BBA Ammon (\APACyear1999) \APACinsertmetastarLigorria1999{APACrefauthors}Ligorría, J\BPBIP.\BCBT \BBA Ammon, C\BPBIJ. \APACrefYearMonthDay1999. \BBOQ\APACrefatitleIterative Deconvolution and Receiver-Function Estimation Iterative deconvolution and receiver-function estimation.\BBCQ \APACjournalVolNumPagesBulletin of the Seismological Society of America8951395–1400. {APACrefDOI} 10.1785/BSSA0890051395 \PrintBackRefs\CurrentBib
Moorkamp \BOthers. (\APACyear2007) \APACinsertmetastarMoorkamp2007{APACrefauthors}Moorkamp, M., Jones, A\BPBIG.\BCBL \BBA Eaton, D\BPBIW. \APACrefYearMonthDay2007. \BBOQ\APACrefatitleJoint Inversion of Teleseismic Receiver Functions and Magnetotelluric Data Using a Genetic Algorithm: Are Seismic Velocities and Electrical Conductivities Compatible? Joint inversion of teleseismic receiver functions and magnetotelluric data using a genetic algorithm: Are seismic velocities and electrical conductivities compatible?\BBCQ \APACjournalVolNumPagesGeophysical Research Letters3416. {APACrefDOI} 10.1029/2007GL030519 \PrintBackRefs\CurrentBib
Muir \BBA Tsai (\APACyear2019) \APACinsertmetastarMuir2019{APACrefauthors}Muir, J\BPBIB.\BCBT \BBA Tsai, V\BPBIC. \APACrefYearMonthDay2019. \BBOQ\APACrefatitleGeometric and Level Set Tomography Using Ensemble Kalman Inversion Geometric and level set tomography using ensemble kalman inversion.\BBCQ \APACjournalVolNumPagesGeophysical Journal International. {APACrefDOI} 10.1093/gji/ggz472 \PrintBackRefs\CurrentBib
Plett (\APACyear2006) \APACinsertmetastarPlett2006{APACrefauthors}Plett, G\BPBIL. \APACrefYearMonthDay2006. \BBOQ\APACrefatitleSigma-Point Kalman Filtering for Battery Management Systems of LiPB-Based HEV Battery Packs: Part 2: Simultaneous State and Parameter Estimation Sigma-point kalman filtering for battery management systems of lipb-based hev battery packs: Part 2: Simultaneous state and parameter estimation.\BBCQ \APACjournalVolNumPagesJournal of Power Sources16121369–1384. {APACrefDOI} 10.1016/j.jpowsour.2006.06.004 \PrintBackRefs\CurrentBib
Sambridge (\APACyear1999) \APACinsertmetastarSambridge1999{APACrefauthors}Sambridge, M. \APACrefYearMonthDay1999. \BBOQ\APACrefatitleGeophysical Inversion with a Neighbourhood Algorithm—II. Appraising the Ensemble Geophysical inversion with a neighbourhood algorithm—II. Appraising the ensemble.\BBCQ \APACjournalVolNumPagesGeophysical Journal International1383727–746. {APACrefDOI} 10.1046/j.1365-246x.1999.00900.x \PrintBackRefs\CurrentBib
Sambridge \BOthers. (\APACyear2013) \APACinsertmetastarSambridge2013{APACrefauthors}Sambridge, M., Bodin, T., Gallagher, K.\BCBL \BBA Tkalčić, H. \APACrefYearMonthDay2013. \BBOQ\APACrefatitleTransdimensional Inference in the Geosciences Transdimensional inference in the geosciences.\BBCQ \APACjournalVolNumPagesPhilosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences371198420110547. {APACrefDOI} 10.1098/rsta.2011.0547 \PrintBackRefs\CurrentBib
Stuart (\APACyear2010) \APACinsertmetastarStuart2010{APACrefauthors}Stuart, A\BPBIM. \APACrefYearMonthDay2010. \BBOQ\APACrefatitleInverse Problems: A Bayesian Perspective Inverse problems: A bayesian perspective.\BBCQ \APACjournalVolNumPagesActa Numerica19451–559. {APACrefDOI} 10.1017/S0962492910000061 \PrintBackRefs\CurrentBib
Tarantola (\APACyear2005) \APACinsertmetastarTarantola2005{APACrefauthors}Tarantola, A. \APACrefYear2005. \APACrefbtitleInverse Problem Theory and Methods for Model Parameter Estimation Inverse problem theory and methods for model parameter estimation. \APACaddressPublisherSociety for Industrial and Applied Mathematics. {APACrefURL} http://epubs.siam.org/doi/book/10.1137/1.9780898717921 {APACrefDOI} 10.1137/1.9780898717921 \PrintBackRefs\CurrentBib
Wang \BOthers. (\APACyear2019) \APACinsertmetastarWang2019{APACrefauthors}Wang, J., Yang, D., Jing, H.\BCBL \BBA Wu, H. \APACrefYearMonthDay2019. \BBOQ\APACrefatitleFull Waveform Inversion Based on the Ensemble Kalman Filter Method Using Uniform Sampling without Replacement Full waveform inversion based on the ensemble kalman filter method using uniform sampling without replacement.\BBCQ \APACjournalVolNumPagesScience Bulletin645321–330. {APACrefDOI} 10.1016/j.scib.2019.01.021 \PrintBackRefs\CurrentBib