Statistical estimation

Let $\bold{d}$ be the vector of observed data, and $\bold{m}$ be the ideal underlying model. The regularized data represent the model estimate $<\!\!\bold{m}\!\!\gt$. Taking into account the lack of information about $\bold{m}$, we can treat both $\bold{m}$ and $\bold{d}$ as random vectors and approach the problem of finding $<\!\!\bold{m}\!\!\gt$ statistically.

For any two random vectors $\bold{x}$ and $\bold{y}$, let us denote by $\bold{C}_{xy}$ the mathematical expectation of the random matrix $\bold{x}\,\bold{y}^T$, where $\bold{y}^T$ denotes the adjoint of y. Analogously, $\bold{C}_{x}$ will denote the mathematical expectation of $\bold{x}\,\bold{x}^T$. For zero-mean vectors, the matrices $\bold{C}_{x}$ and $\bold{C}_{xy}$ correspond to covariances. In a more general case, they are second-moment statistics of the corresponding random processes.

Applying the Gauss-Markoff theorem, one can obtain an explicit form of the estimate $<\!\!\bold{m}\!\!\gt$ under three very general assumptions Liebelt (1967):

1.

The estimate has a linear relationship with the input data:  

$$<\!\!\bold{m}\!\!\gt = \bold{A}\,\bold{d}\;,$$ (1)

where $\bold{A}$ is a linear operator.

2.

The estimate corresponds to the minimum of $\bold{C}_e = E\left(\bold{e}\,\bold{e}^T\right)$, where E is the mathematical expectation and $\bold{e}$ denotes the model error $\bold{e} = <\!\!\bold{m}\!\!\gt - \bold{m}$. For unbiased estimates (zero mathematical expectation of $\bold{e}$), the matrix $\bold{C}_e$ corresponds to the model error covariance. Although we do not make any explicit assumptions about the statistical distribution of the error, minimizing $\bold{C}_e$ is particularly meaningful in case of normal (Gaussian) distributions Tarantola (1987).

3.

The square matrix $\bold{C}_d$ is invertible.

Doing a simple algebraic transformation, we find that

$$\begin{eqnarray} \nonumber \bold{C}_e & = & E\left[\left(<\!\!\bold{m}\!\!\gt -... ...T - \bold{C}_{md}\,\bold{C}_d^{-1}\,\bold{C}_{md} + \bold{C}_m\;.\end{eqnarray}$$ (2)

It is evident from equation () that C_e will be minimized when $\bold{A} = \bold{C}_{md}\,\bold{C}_d^{-1}$. This leads immediately to the Gauss-Markoff result  

$$<\!\!\bold{m}\!\!\gt = \bold{C}_{md}\,\bold{C}_d^{-1}\,\bold{d}\;.$$ (3)

Equation () has fundamental importance in different data regularization schemes. With some slight modifications, it appears as the basis for such methods as optimal interpolation in atmospheric data analysis Daley (1991); Gandin (1965), least-squares collocation in geodesy Moritz (1980), and linear kriging in petroleum and mining engineering Hohn (1999); Journel and Huijbregts (1978). In order to apply formula () in practice, one needs first to get an estimate of the matrices $\bold{C}_{md}$ and $\bold{C}_d$. In geostatistics, the covariance matrices are usually chosen from simple variogram models Deutsch and Journel (1997).

Unfortunately, a straightforward application of the Gauss-Markoff formula () is computationally unaffordable for typical seismic data applications. If the data vector contains N parameters, a straightforward application will lead to an N by N matrix inversion, which requires storage proportional to N² and a number of operations proportional to N³. Although the data can be divided into local patches to reduce the computational requirements for an individual patch, the total computational complexity is still too high to be affordable for the values of N typical in 3-D seismic exploration (N as high as 10¹⁰).

We can take two major theoretical steps to reduce the computational complexity of the method. The first step is to approximate the covariance matrices with sparse operators so that the matrix multiplication is reduced from N² operations to something linear in N. The second step is to approach model estimation as an optimization problem and to use an iterative method for solving it. The goal is to obtain a reasonable model estimate after only a small number of iterations.