Constrained Full Waveform Inversion for Borehole Multicomponent Seismic Data

Abstract

Full-waveform inversion for borehole seismic data is an ill-posed problem and constraining the problem is crucial. Constraints can be imposed on the data and model space through covariance matrices. Usually, they are set to a diagonal matrix. For the data space, signal polarization information can be used to evaluate the data uncertainties. The inversion forces the synthetic data to fit the polarization of observed data. A synthetic inversion for a 2D-2C data estimating a 1D elastic model shows a clear improvement, especially at the level of the receivers. For the model space, horizontal and vertical spatial correlations using a Laplace distribution can be used to fill the model space covariance matrix. This approach reduces the degree of freedom of the inverse problem, which can be quantitatively evaluated. Strong horizontal spatial correlation distances favor a tabular geological model whenever it does not contradict the data. The relaxation of the spatial correlation distances from large to small during the iterative inversion process allows the recovery of geological objects of the same size, which regularizes the inverse problem. Synthetic constrained and unconstrained inversions for 2D-2C crosswell data show the clear improvement of the inversion results when constraints are used.

1. Introduction

Full-waveform inversion of seismic data allows one to obtain an image of the subsurface through the determination of a certain number of physical parameters. Borehole data (VSP and crosswell) are of special interest because of the geometry of acquisition, which provides more informative signals on the medium properties than surface seismic data (high frequency signal, energetic P-S conversions). However, the lack of seismic data redundancy (e.g., few shots) renders the inversion problem underdetermined. In general, full-waveform inversion is an ill-posed problem, in the sense that an infinite number of models match the data [1]. It is classically solved with local optimization schemes and is therefore strongly dependent on the starting model definition. This starting model should predict arrival times with errors less than half of the period to cancel the cycle-skipping ambiguity [1]. The multiscale strategy performed by moving from low to high frequencies during the inversion allows reduction of the nonlinearities and cycle skipping issues of the inversion and helps convergence toward the global minimum. Regularizations are conventionally applied to the inversion in the model space to make it better posed [2,3]. Tikhonov and Arsenin [4] have proposed a regularization strategy within the optimization step to find the smoothest model that explains the data. Preconditioning techniques acting as a smooth operator on the model update [5] may add strong prior features of the expected structure through directive Laplace preconditioning, as in Guitton et al. [6]. Regularization schemes that preserve edges and contrasts have also been developed for specific full waveform inversion applications through a ${\ell }_1$ model penalty [7] or through a multiplicative regularization [8]. Regularization can also be expressed in the curvelet or wavelet domains [9,10]. In such domains, the ${\ell }_1$ norm minimization is generally preferred for the model term penalty because it ensures sparsity in the model space. All the previous regularization techniques allow stabilization of the inversion scheme by assuming a particular representation or structure of the velocity model (smoothness, sparsity, and so on). However, geological information, taking into account prior model information, is generally not used in classical full waveform inversion implementation; however, there are some examples [11,12]. In addition, most of the weighting operators or covariance matrices associated with the model parameters are set to identity or in the best cases to a diagonal matrix. The same applies to the weighting operators or covariance matrices associated with the data.

Weighting operators on the data and models or inverse of the data and model covariance operators for least squares in the frame of the Bayesian formulation [3] are introduced in the misfit function. In this paper, we will show the benefit of introducing covariance matrices in the data and model space. On the one hand, we will illustrate the impact on constraining the seismic inverse problem on the data space for the case of a synthetic two-component borehole seismic data by using the polarization analysis to fill a block diagonal data space covariance matrix. On the other hand, we will illustrate the benefit of constraining the seismic inverse problem on the model space by performing a crosswell synthetic experiment. The horizontal and vertical spatial correlations using Laplace distribution are used to fill the model space covariance matrix in order to introduce an a priori solution, favoring a tabular medium whenever it doesn’t contradict the information contained in the data.

2. Constrained Least Squares Inversion

Using the probabilistic formalism developed by Tarantola [3] and considering only the data space term in the equation, we rewrite the misfit function in the vicinity of the current model $m_n$:

$$\begin{array}{ccc}S\left(m_n+\delta m\right)& =& {\left(\delta d_n-G_n\delta m\right)}^T{C}_D{}^{-1}\left(\delta d_n-G_n\delta m\right)\\ & & +{\left(\Delta m_n-\delta m\right)}^T{C}_M{}^{-1}\left(\Delta m_n-\delta m\right),\end{array}$$

(1)

where:

• $\delta m$ is the perturbation in the vicinity of the current model $m_n$,
• $\delta d_n$ are the data residuals for the model $m_n$,
• $\Delta m_n$ is the difference between the current model $m_n$ and the a priori model $m_{prior}$
• $G_n$ is the linear function tangent to $g$ at the model $m_n$,
• $g$ is the function mapping the model space $m\in \mathit{M}$ into the data space $d=g\left(m\right)\in \mathit{D}$
• $C_D$ is the covariance matrix on the data space.
• $C_M$ is the covariance matrix on the model space (defining the Gaussian a priori probability density).

Applied to the waveform inversion problem of estimating the elastic parameters and the density of the earth, the minimization of the misfit function can be solved by iterative gradient methods.

The iterative process of a nonlinear inversion can be summarized as follows [3]: for a given iteration $k$, from the current model $m_k$ (i.e., the discretized physical fields characterizing the medium), we perform simulations of the wave in order to obtain the synthetic data $d_k^{cal}=g\left(m_k\right)$. Thus, comparing them to observed data $d^{obs}$, we obtain the residuals $\delta d_k=d^{obs}- d_k^{cal}$. These residuals are weighted ${\widehat{\delta d}}_k=C_D^{-1}\delta d_k$ (the hat denotes the dual space) and back-propagated using the wave propagation equation again to obtain the gradient in the dual model space ${\widehat{\mathsf{\gamma }}}_k= G_k^T{C}_D^{-1}\delta d_k$ (where $G_k$ is the derivative of the function $g$ over the model space at the point $m_k$ and where $T$ denotes the transpose operation) and consequently, the gradient in the model space ${\mathsf{\gamma }}_k=C_M{\widehat{\mathsf{\gamma }}}_k$. The gradient indicates the steepest ascent direction of the misfit function in the model space. We modify this direction using the conjugate gradient algorithm ${\mathsf{\phi }}_k={\mathsf{\gamma }}_k+ {\alpha }_k{\mathsf{\phi }}_{k-1}$, where ${\mathsf{\phi }}_k$ is the conjugate gradient and ${\alpha }_k$ is defined from the gradient [13]. Finally, using an additional simulation, we optimize the step length ${\mu }_k$ [14] in the upgrading equation $m_{k+1}=m_k-{\mu }_k{\mathsf{\phi }}_k$ and we update the model.

3. Constraints on the Data Space

The quantification of uncertainties in seismic data is often neglected, because in many cases it is a difficult task. For instance, data errors generated by receiver response are correlated in time, which means accounting for them is not straightforward. To illustrate the possible contribution of constraints provided by the analysis of uncertainties of data, we have performed two inversions from the synthetic data experiment. The first inversion is done without any data constraints, whereas in the second one we have incorporated polarization wave analysis in the covariance matrix $C_D$. This matrix treats the uncertainties on particle velocities according to the eigenbasis of the particle motion.

For the case of multi-component data, the uncertainty of a component may be correlated to the uncertainty of the other components. For each time for each of the receiver components, we can define a 2 × 2 matrix for 2-component data (as in the example bellow), or a 3 × 3 matrix for 3-component data. In that case, the covariance matrix on the data space $C_D$ is not a diagonal matrix. In order to take into account polarization, we consider that the uncertainties are proportional to the local cross-correlation matrix of the signal (the data) and that the signal is locally stationary (over 2 periods). The cross-correlation matrix between components $i$ and $j$ of the signal $s_{r,t}$ is computed with zero lag ($\tau =0$) for any receiver $r$ and any time $t$:

$$C_{r,t}^{ij}\left(\tau \right)={{\displaystyle \int }}_{-\infty }^{+\infty }{s}_{r,t}^i\left(u\right)s_{r,t}^j\left(u+\tau \right)du,$$

(2)

where:

$$s_{r,t}^i\left(u\right)=d_r^i\left(u\right) w_t\left(u\right)$$

with $d_r^i\left(u\right)$ being the trace for the component $i$ of the receiver $r$ of the observed data and $w_t\left(u\right)$ a Hamming time window centered on time $t$ with a length of typically 2 periods. One can use other time tapering windows classically used in spectral analysis (Blackman, Nuttall, etc).

3.1. Cross-Correlations and Polarization

The polarization of the multicomponent signal is estimated from this matrix $C_{r,t}^{ij}\left(\tau =0\right)$ and the matrix of the uncertainties has the same eigenvectors (polarization directions) than the cross-correlation matrix. The uncertainty has the same shape as the polarization of the signal. The ratios among the eigenvalues are the same, in other words, the rectilinearity and planarity of the polarization [15] is conserved.

The polarized signal for rectilinear polarization can be modeled by:

$$d\left(t\right)=a·s\left(t\right)+n\left(t\right),$$

(3)

where $d\left(t\right)$ is the 3C signal, $a$ is the polarization vector, $s\left(t\right)$ the scalar signal along the polarization vector, and $n\left(t\right)$ the 3C Gaussian noise.

In the frequency domain, the equation can be expressed as:

$$\widehat{d}\left(f\right)=a·\widehat{s}\left(f\right)+\widehat{n}\left(f\right),$$

(4)

where the hat denotes the frequency domain.

Let us consider the spectral matrix of the 3C signal defined by:

$$\begin{array}{lll}\widehat{S}\left(f\right)& =& \overline{\widehat{d}}\left(f\right)·{\widehat{d}}^T\left(f\right)\\ & =& \overline{\left(a\widehat{s}+\widehat{n}\right)}·{\left(a\widehat{s}+\widehat{n}\right)}^T\\ & =& a·a^T·\overline{\widehat{s}}\widehat{s}+\overline{\widehat{\mathit{n}}}·{\widehat{n}}^T\\ & =& a·a^T{\left|\widehat{s}\left(f\right)\right|}^2+\widehat{N}\left(f\right),\end{array}$$

(5)

where the upper bar denotes the complex conjugate, $T$ denotes the transpose operator and $\widehat{N}\left(f\right)$ being the spectral matrix of the noise. Note that $\overline{\left(a\widehat{s}\right)}·{\widehat{n}}^T=\overline{\widehat{n}}·a^T\widehat{s}=0$ because noise and signal are considered uncorrelated. The polarization vector is an eigenvector of the spectral matrix (for isotropic noise) for any frequency $f$:

$$\begin{array}{lll}\widehat{S}·a& =& (a·a^T{\left|\widehat{s}\right|}^2+\widehat{N})a\\ & =& {\left|\widehat{s}\right|}^{2}a·a^{T}a+\widehat{N}a\\ & =& {\left|\widehat{s}\right|}^{2}a+({\left|\widehat{n}\right|}^{2}I)a\\ & =& ({\left|\widehat{s}\right|}^2+{\left|\widehat{n}\right|}^2)a.\end{array}$$

(6)

When integrating the spectral matrix over $ℝ$, one obtains the cross-correlation matrix $C\left(\tau =0\right)$ with zero lag with the same eigenvectors.

The structure of the covariance matrix on the data space $C_{\mathrm{D}}$ is derived from the previous equations. For instance, for a 2C receiver, a given source or receiver couple, and at any time $t$, the corresponding 2 × 2 submatrix is defined as:

$$C_{\mathrm{D}}=\mathit{V}·{\mathsf{\sigma }}_{\mathrm{D}}{}^2\left(\begin{array}{cc}1& 0\\ 0& \frac{{\mathsf{\lambda }}_{\mathrm{II}}}{{\mathsf{\lambda }}_{\mathrm{I}}}\end{array}\right)·{\mathit{V}}^{-1}$$

where $\mathit{V}$ is the eigenvectors matrix (column), ${\mathsf{\sigma }}_{\mathrm{D}}$ is the standard deviation of the (isotropic) noise for this couple, and ${\mathsf{\lambda }}_{\mathrm{I}}$ and ${\mathsf{\lambda }}_{\mathrm{II}}$ are the max and min of the two positive eigenvalues of the spectral matrix, respectively.

Let us consider the sub-matrix for the first data 2C-trace of the covariance matrix on the data space:

$$C_{\mathrm{D}}\left(rec=1\right)\propto C_{\mathrm{D}}\left(r=1,t,\tau =0\right)$$

$$\left[\begin{array}{cccccccc}{\sigma }_{1,X}^2& {\rho }_1{\sigma }_{1,X}{\sigma }_{1,Z}0& 0& 0& \cdots & \cdots & 0& 0\\ {\rho }_1{\sigma }_{1,X}{\sigma }_{1,Z}& {\sigma }_{1,Z}^2& 0& 0& \cdots & \cdots & 0& 0\\ 0& 0& {\sigma }_{2,X}^2& {\rho }_2{\sigma }_{2,X}{\sigma }_{2,Z}& \cdots & \cdots & 0& 0\\ 0& 0& {\rho }_2{\sigma }_{2,X}{\sigma }_{2,Z}& {\sigma }_{2,Z}^2& \cdots & \cdots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & \ddots & ⋮& ⋮\\ ⋮& ⋮& ⋮& ⋮& \ddots & \ddots & ⋮& ⋮\\ 0& 0& 0& 0& \cdots & \cdots & {\sigma }_{nt,Z}^2& {\rho }_{nt}{\sigma }_{nt,X}{\sigma }_{nt,Z}\\ 0& 0& 0& 0& \cdots & \cdots & {\rho }_{nt}{\sigma }_{nt,X}{\sigma }_{nt,Z}& {\sigma }_{nt,Z}^2\end{array}\right]$$

where ${\rho }_i$, ${\sigma }_{i,X}$, and ${\sigma }_{i,Z}$ are, respectively, the correlation coefficient, the standard deviation for the $X$ component, and the standard deviation for the $Z$ component for a given time sample $i$, and $nt$ being the total number of time samples. The matrix $C_{\mathrm{D}}$ is a block diagonal matrix and its inverse (i.e., $C_{\mathrm{D}}^{-1}$, is another block diagonal matrix), composed of the inverse of each block.

3.2. An Offset VSP Synthetic Example

In this numerical example, the acquisition geometry is an offset VSP with a unique source at 500 m offset and the antenna is between 580 and 970 m depth in a vertical well (see Figure 1 for more details). The isotropic elastic model depends only on the depth (1D model). The “true” model is defined in Figure 1. The horizontal and vertical components of the “observed” data obtained from the true model are displayed in Figure 1. The direct problem (seismic modeling) is based on discretizing the wave equation by the finite differences method (FDM) [16].

The inverted parameters are the P-wave velocity, the S-wave velocity, and the density. The starting models are simple, the water layer is considered as known, and the other layers are replaced by a constant vertical gradient for each parameter (gray solid lines in Figure 2). The test consists of comparing the results of inversion when data uncertainties are isotropic and when the polarization is included in the data uncertainties.

As shown in Figure 2, the inversion results exhibits in both cases a good fit in the receiver zone, while the estimation is not good in the upper part of the model. One can clearly note that the S-wave velocity model is recovered with more detail than the P-wave velocity model. This is due to the fact that converted S-waves are present in borehole seismic data and also because the S-waves have a shorter wavelength than P-waves, hence, a better spatial resolution. The final residuals for both inversion experiments are small; a few percent of the misfit. The inversion results when polarization is included are significantly better than for the isotropic uncertainties. More precisely, for the case of polarization, both P-wave and S-wave estimated velocity fields are well fitted in the antenna zone.

4. Constraints on the Model Space

We can distinguish two kinds of geological a priori information that can be incorporated in an inversion process. On the one hand, information obtained by measurements (well logging), that we consider as “objective” information, allows the definition of an a priori model at the well vicinity. On the other hand, information coming from geological interpretations, such as geological layer dips, that we qualify as “subjective” information, can be simplified and introduced into the covariance matrix $C_M$.

4.1. Quantification of Number of Degrees of Freedom

In order to evaluate the importance of constraints on an inversion, it is necessary to quantify the number of degrees of freedom. This number can be estimated from the set of parameters as a function of existing correlations among them and independently of their standard deviations. The number $N$ of free parameters is equal to the number $n$ of parameters in the absence of correlations, and it is reduced to 1 in the case of perfect correlations.

4.1.1. General Case

Let us consider a 1D Gaussian random field $X$, the domain is discrete and finite consisting of a set of $n$ points; this field is equivalent to the set of random variables $\left(X_1, X_2,\cdots ,X_n\right)$. It can be characterized by the following probability density function (or p.d.f.):

$$f\left(x\right)=\frac{1}{{\left(2\pi \right)}^{n/2}\det {\left(C\right)}^{1/2}}\exp \left(-\frac{1}{2}{\left(x-m\right)}^T{C}^{-1}\left(x-m\right)\right),$$

(7)

where $x$ is a realization, $m$ is the expectation of the field $\left(\mathrm{E}\left(X\right)=m\right)$, and $C$ is the covariance matrix of the field.

The covariance matrix can be rewritten with the following form $=SRS$, where $R$ is the correlation matrix and $S$ the diagonal matrix of the standard deviations.

If there is no correlation among the random variables $X_i$, then $R=I$ and the covariance matrix is diagonal and the number of free parameters is then $n$. If a correlation exists, we then use the $LU$ decomposition method to estimate the number of free parameters (see Appendix A and Appendix B for more details). We can write $R=M{M}^T$ where $M$ is the left lower triangle matrix of the $LU$ decomposition of the matrix $R$. The number of free parameters is then reduced to:

$$N=trace\left(M\right).$$

(8)

4.1.2. Spatial Correlation on the Model Space

For illustration purposes, let us consider the special case where the spatial correlations are described by an exponential model. Thus, for any pair of points $M_1$ and $M_2$ (belonging to the same horizontal line), the correlation $\rho$ between the pair of random variables $X_1$ and $X_2$ corresponding to the points $M_i$ can be expressed as a function of the distance:

$$\rho \left(X_1,X_2\right)=\exp \left(-\frac{d\left(M_1,M_2\right)}{r}\right),$$

(9)

where $d(·)$ is a distance and $r$ is the range of the correlation.

By setting the variance for each variable to be constant ${\sigma }_i={\sigma }_1$, the p.d.f. of Equation (7) reduces to,

$$f\left(x\right)=\frac{1}{{\left(2\pi \right)}^{n/2}·{\sigma }^n\sqrt{\det \left(R\right)}}\exp \left(-\frac{1}{2}{\left(\frac{x-m}{\sigma }\right)}^T{R}^{-1}\left(\frac{x-m}{\sigma }\right)\right),$$

(10)

This simplifies, as we will show below, the study of the relation between the matrix $R$ and the degrees of freedom using an exponential correlation function.

Let us consider the same Gaussian field $X$, but this time defined on a defined $1\mathrm{D}$ domain regularly sampled (as for a grid) containing $n$ points. The distances between points are multiples of the step $\Delta x$ and the correlation matrix has the following form (for ordered points):

$$R=\left(\begin{array}{ccccccc}1& a& a^2& a^3& a^4& \cdots & a^n\\ a& 1& a& a^2& a^3& \cdots & a^{n-1}\\ ⋮& ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ a^n& a^{n-1}& a^{n-2}& a^{n-3}& a^{n-4}& \cdots & 1\end{array}\right),$$

(11)

where

$$a=\exp \left(-\frac{\Delta x}{r}\right)\le 1.$$

(12)

One interesting property of the $R$ matrix is that the $M$ matrix of the $LU$ decomposition has a simple form

$$M=\left(\begin{array}{ccccccc}1& 0& 0& 0& 0& \cdots & 0\\ a& b& 0& 0& 0& \cdots & 0\\ a^2& ab& b& 0& 0& \cdots & 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ a^n& a^{n-1}b& a^{n-2}b& a^{n-3}b& a^{n-4}b& \cdots & b\end{array}\right),$$

(13)

with $b=\sqrt{1-a^2}$. This means that, in terms of a sequential simulation, the $n-1$ variables, $X_2$ to $X_n$, have the same degree of freedom: $b$. The total degrees of freedom $N$ of the field $X$ is then:

$$N=1+\left(n-1\right)·b,$$

(14)

where for $n$ parameters (or $n$ grid points), the number of free parameters, considered as the degrees of freedom, of the inversion is $N$. Figure 3 shows the relation between the range of the correlation (expressed in number of steps, $\Delta x$) and the values of $b$.

4.2. Seismic Crosswell Numerical Experiment Example

This numerical example is an illustration how to incorporate the “subjective” geological information that our medium is mainly tabular but might contain complex geological structures using the covariance matrix on the model space.

4.2.1. Description of the Crosswell Experiment and Observed Data

The reference subsurface models are shown in Figure 4. The word “reference” indicates that this model is used to generate the synthetic seismic data used as observed data (see Figure 5) during the inversion process. It is also the model to be retrieved by applying the full-wave inversion method to the observed data. The medium for wave propagation is elastic isotropic. The parameters to be inverted are discretized physical fields: vertical P-wave and S-wave velocities. For simplicity, the density ($\rho$) is not inverted—the field is constant with $\rho =2500 \mathrm{kg}·{\mathrm{m}}^{-3}$ and it fixed at the correct value in the different inversions.

4.2.2. Multiscale Constrained Inversion

The inversion process consists of four successive inversions using the same observed data but with decreasing correlation ranges in horizontal and vertical directions; therefore, it is a multiscale inversion and each inversion is for a given “scale”. The model parameter results of each inversion are used as initial model parameters for the next inversion.

Uncertainties on the model space are required for the definition of $C_M{}^{-1}$. It can be defined by specifying the spatial correlation ranges and the uncertainties on the parameters. Here the uncertainties on the parameters are kept constant: ${\sigma }_{Vp}=120 \mathrm{m}/\mathrm{s}$ and ${\sigma }_{Vs}=100 \mathrm{m}/\mathrm{s}$. The ranges of the horizontal and vertical correlations do not depend on the parameters but on the spatial location. The horizontal correlation depends on depth and both horizontal and vertical ranges depend on the inversion scale. We define three types of a priori regions with respect to depth (see

Table 1. A priori geological region type with respect to depth.

Depth (m)	Region Type	Number of Vertical Points
800–1120	Quasi 1D	161
1120–1180	Transition	29
1180–1340	2D	81
1340–1400	Transition	29
1400–1600	Quasi 1D	101

): “Quasi 1D” with large correlation ranges, “2D” with short correlation ranges (where we expect to have complex geological objects, see Figure 4), and “Transition”, with correlation ranges linearly varying with depth to match the ranges of the two adjacent regions.

Table 2. Horizontal and vertical correlation ranges for the different a priori region types and inversion scales.

Region Type	Correlation Ranges (m)	Inversion Scale (b)	Inversion Scale (c)	Inversion Scale (d)	Inversion Scale (e)
Quasi 1D	$(r_x, r_z)$	(1000, 24)	(1000, 10)	(1000, 6)	(1000, 2)
2D	$(r_x, r_z)$	(320, 24)	(80, 10)	(20, 6)	(5, 2)

lists the different correlation ranges associated with the four successive inversion scales and for the a priori region types. The four successive inversion scales are denoted by (b), (c), (d), and (e) in relation to the results displayed in Figure 6 and Figure 7. The correlation ranges are gradually decreasing with the successive inversion scales in order to solve first the large wavelengths in the model before resolving the small wavelengths (i.e., the small model structures).

The results of the multiscale inversion process for the P-wave velocity field are displayed in Figure 6. The P-wave velocity starting model is obtained by a low-frequency full waveform inversion. It is also used as the starting model for the unconstrained inversion; such starting models can also be obtained from traveltime inversion, especially for crosswell data. The estimated P-velocity field is obtained from the first inversion scale, i.e., using correlation ranges of 320 m and 24 m in the 2D region, respectively, for the horizontal and vertical correlations. One can notice that the 2D structure cannot be resolved, as the horizontal correlation range is larger than the size of the triangular object. The 1D structures are partially retrieved. For the inversion scale (c), the estimated P-velocity field is, obtained from the inversion scale (b), i.e., using correlation ranges of 80 m and 10 m in the 2D region, respectively, for the horizontal and vertical correlations. One can notice that the 2D structure start to be resolved as both correlation ranges are smaller than the size of the triangle. The 1D structures are well estimated, with some vertical smoothening due to the vertical correlation range of 10 m. For the inversion scale (d), the estimated P-velocity field is obtained from the third inversion scale, i.e., using correlation ranges of 20 m and 6 m in the 2D region, respectively, for the horizontal and vertical correlations. One can notice that the 2D structure is getting better spatial resolution, for both the triangle and the lower zigzag interface. The 1D structures are well resolved. For the inversion scale (e), the estimated P-velocity field is obtained from the fourth and last inversion scale, i.e., using correlation ranges of 5 m and 2 m in the 2D region, respectively, for the horizontal and vertical correlations. One can notice that the delineation of the 2D structure is close to the reference model. As expected, the estimated field outside the zone between the two wells is not reliable.

The results of the multiscale inversion process for the S-wave velocity field are displayed in Figure 7. The results are similar to the P-wave velocity field (Figure 6). The main difference is clear at the fourth scale (e): the spatial resolution is almost perfect for the S-wave due to smaller wavelength content associated to this mode. However, the inversion process provides an overestimation of velocity contrast near the interfaces and some small oscillations in the estimated field clearly visible for the light blue domain above the central triangle. It means that the number of degrees of freedom is too high regarding the amount of information provided by the data.

The misfit function for the successive multiscale inversion continuously decreases from one inversion scale to another (see

Table 3. The misfit and iteration number for the different inversion scales.

Inversion Parameters	Inversion Scale (b)	Inversion Scale (c)	Inversion Scale (d)	Inversion Scale (e)
Initial Misfit	45.5%	27.9%	6.55%	0.81%
Final Misfit	27.9%	6.55%	0.81%	0.14%
Iteration #	30	30	50	50

).

The number of degrees of freedom per point are provided in

Table 4. Number of degrees of freedom (DoF) per point for the successive multiscale inversions.

Region Type	Inversion Scale (b)	Inversion Scale (c)	Inversion Scale (d)	Inversion Scale (e)
Quasi 1D	1919	2812	3416	4554
Transition	267	468	622	857
2D	1022	2922	6792	15,725
All	3208	6202	10,830	21,136
DoF/point	2.8%	5.5%	9.6%	18.8%

for the different inversion scales. These numbers of freedoms are calculated from Equation (14). The ratio of the number of degrees with respect to the number of model parameters (i.e., grid points), in other words the number of degrees of freedom per point, is provided. It is increasing slower for the quasi 1D region than for the 2D region because only the vertical correlation range is decreasing, while both horizontal and vertical ranges are decreasing in the 2D region. These ratios are pointing out that the degrees of freedom are significantly lower in the early scales of the multiscale inversion, allowing better constraint of the problem and improving the stability by defining an overdetermined system. During the last scale stages, the constraints are relaxed, allowing one to obtain a better resolution while avoiding artifacts in the estimated fields. Even for the last scale, the ratio indicates about 5 times less free parameters than in the unconstrained inversion.

For comparison purposes with the results of the constrained inversion, we performed an unconstrained inversion. In the Figure 8, we display the results of the multiscale constrained inversion with respect to the unconstrained inversion for the same misfit of 0.3%. We can notice that the unconstrained inversion results are noisier than the constrained inversion. The main artifacts can be found near the sources (classical problem in borehole seismic full waveform inversion), below the triangle structure, above the source or receiver zones (migration smile type artifacts), and in the constant gradient 1D layers. These artifacts in the estimated fields are the consequence of the instability of the least-squares inversion problem when the problem is not correctly determined for all parameters, in other words, when it is partially underdetermined. The number of degrees of freedom is too large regarding the information provided by the data. Reducing the number of degrees of freedom of the problem taking into account prior geological information (i.e., constraining the fields using spatial statistic where we have more prior information) allows one to overcome this problem, reducing drastically the number and the magnitude of the artifacts.

5. Conclusions

In this paper, we have illustrated how to introduce, in the inverse problem, constraints consistent with the least squares formalism (through covariance matrices) on both the data and the model space.

Constraining the seismic inverse problem in the data space is not straightforward but when possible, it is a very powerful tool. In the case of a multicomponent borehole seismic data, using the polarization in the inversion process allowed better recovery of the sharpness of the interfaces at the level of the receivers.

Constraining the seismic inverse problem in the model space can be achieved by defining a priori model parameters and evaluating the model uncertainties associated with this model. Horizontal and vertical spatial correlation using Laplace distribution can be used to fill the model space covariance matrix in order to introduce an a priori information favoring the tabular region versus more complex regions, by varying accordingly the ranges of these correlations. The inversion favors solutions are consistent with our a priori, whenever it does not contradict the information contained in the data. Moreover, by adopting a multiscale type of inversion by relaxing the correlation ranges, it regularizes the inverse problem.

The incorporation of all our a priori knowledge of the parameters and all statistical studies on the data, allows not only the algorithmic stabilization of the inversion process, but also the reduction of the solution set for an underdetermined problem, with the purpose being not necessarily to converge quickly towards a good model (in term of residuals), but to prospect regions of the model space populated by models that are sensible, a priori, and also yielding the lowest possible misfit.