5. Simulation of process trajectories

The classical method for simulating a stationary, centered Gaussian stochastic process is based on the spectral representation of the processes, cf. for example [bib11]. The spectral representation as well as the numerical simulation algorithm are described below.

We treat the general case of a Gaussian vector process \(Y(t)\in {ℝ}^{M},t\in [\mathrm{0,}T]\) characterized by its DSP matrix \({S}_{Y}(\omega )\in {\mathrm{Mat}}_{ℂ}(M,M)\). This method can be extended to the case of non-stationary processes to DSP scalable \({S}_{X}(\omega ,t)\in {\mathit{Mat}}_{ℂ}(M,M)\) if the process is slowly changing.

5.1. Discretization

Numerical simulation requires time domain discretization. We use a time discretization with a constant \(\Delta t\) step. The cutoff frequency is deduced directly from this value: \({\Omega }_{C}=\pi /(\Delta t)\mathrm{.}\) For discretization by N time steps, the temporal sampling points are denoted \({t}_{j}=j\Delta t,j=\mathrm{0,}\mathrm{...}N-1\). Moreover, the frequency step is defined by \(\Delta \omega =2{\Omega }_{C}/N\) and the duration of the time signal is \(T=\Delta t(N-1)\).

The spectral densities are defined over the interval \(F=[-{\Omega }_{C},+{\Omega }_{C}]\) and the frequency discretization points are chosen as the middle of the tiles:: \({\omega }_{j}=-{\Omega }_{C}+(0.5+j)\Delta \omega ,j=\mathrm{0,}\mathrm{...}N-1\).

5.2. Spectral representation of stationary processes

The spectral representation theorem (cf. [bib11]) says that there is a stochastic vector \(dZ(\omega )\mathrm{\in }{\mathrm{ℂ}}^{M}\) process with orthogonal increments such as:

\(Y(t)\mathrm{=}{\mathrm{\int }}_{{\mathrm{ℝ}}^{M}}{\mathrm{e}}^{i\omega t}dZ(\omega )\) (29)

This process, called the spectral process associated with \(Y\), verifies:

\(E(dZ(\omega )dZ{(\omega \text{'})}^{\ast })=\{\begin{array}{cc}0& \mathit{si}\omega \ne \omega \text{'}\\ {S}_{Y}(\omega )d\omega & \mathit{si}\omega =\omega \text{'}\end{array}\), (30)

where \(E\) is the mathematical expectation operator and \(dZ\) is a process with orthogonal increments.

Since the spectral density matrix \({S}_{Y}(\omega )\) is, by construction, a positive Hermitian matrix for all \(\omega \mathrm{\in }F\), there is a matrix \(L(\omega )\mathrm{\in }{\mathrm{ℂ}}^{M}\) such that:

\({S}_{Y}(\omega )\mathrm{=}L(\omega )L{(\omega )}^{\mathrm{\ast }}\) (31)

Matrix \(L(\omega )\) can be obtained by Cholesky decomposition if \({S}_{Y}(\omega )\) is positive. \(L(\omega )\) is then a lower triangular matrix. In the more general case, decomposition into eigenvectors can be used. If the process is scalar, the problem is reduced to calculating the root of DSP. Process trajectories can then be simulated using the following formula:

\(Y(t)\mathrm{=}\sqrt{\Delta \omega }\mathrm{\Re }e{\mathrm{\sum }}_{j\mathrm{=}0}^{N\mathrm{-}1}L({\omega }_{j}){\chi }_{j}{\mathrm{e}}^{i{\omega }_{j}t}\) (32)

where \({\chi }_{j}\) are complex Gaussian random vectors whose components are independent.

The use of IFFT makes it possible to considerably reduce the numerical cost of simulation. Indeed, we can write the expression (26) in the form

\(Y({t}_{k})\mathrm{=}\sqrt{\Delta \omega }\mathrm{\Re }e({V}_{k}{\mathrm{e}}^{\mathrm{-}i\pi k(1\mathrm{-}1\mathrm{/}N)})\) (33)

where \({V}_{k}\) can be calculated by IFFT:

\({V}_{k}\mathrm{=}{\mathrm{\sum }}_{j\mathrm{=}0}^{N\mathrm{-}1}L({\omega }_{j}){\chi }_{j}{\mathrm{e}}^{2i\pi kj{N}^{\mathrm{-}1}}\)

This algorithm can be used for stationary and non-stationary Gaussian processes with DSP separable. In these cases, the modulated seismic signals are obtained by applying the modulation function so that \(X(t)=q(t)Y(t)\). This is done by the GENE_ACCE_SEISME operator.

5.3. Spectral representation of non-stationary processes at DSP evolving

The spectral representation still applies to non-stationary processes described by their non-separable evolutionary DSP:

\(Y(t)\mathrm{=}{\mathrm{\int }}_{{\mathrm{ℝ}}^{M}}{\mathrm{e}}^{i\omega t}dZ(\omega ,t)\) (34)

provided that DSP evolves slowly over time [bib8,bib6]. There is then a stochastic process \(dZ(\omega ,t)\in {ℂ}^{M}\) with orthogonal increments such as

\(E(dZ(\omega ,t)dZ{(\omega \text{'},t)}^{\mathrm{\ast }})\mathrm{=}\mathrm{\{}\begin{array}{cc}0& \mathit{si}\omega \mathrm{\ne }\omega \text{'}\\ {S}_{Y}(\omega ,t)d\omega & \mathit{si}\omega \mathrm{=}\omega \text{'}\end{array}\) (35)

where \(E\) is the mathematical expectation operator and \(dZ\) is a process with orthogonal increments. Process trajectories can then be simulated using the formula

\(Y(t)\mathrm{=}\sqrt{\Delta \omega }\mathrm{\Re }e{\mathrm{\sum }}_{j\mathrm{=}0}^{N\mathrm{-}1}L({\omega }_{j},t){\chi }_{j}{\mathrm{e}}^{i{\omega }_{j}t}\) (36)

where \({\chi }_{j}\) are complex Gaussian random vectors whose components are independent. The \(L(\omega ,t)\) matrix can still be obtained by Cholesky decomposition if the rank of \({S}_{Y}(\omega ,t)\) is maximum. On the other hand, it is not possible to use the fast Fourier transform as in the stationary case (cf. equation (33)) which increases the calculation times.

The modulated seismic signals are obtained by applying the modulation function so that \(X(t)=q(t)Y(t)\). This is achieved in*Code_Aster* by the GENE_ACCE_SEISME operator.