Velocity model building from seismic reflection data by full‐waveform inversion

Full‐waveform inversion is re‐emerging as a powerful data‐fitting procedure for quantitative seismic imaging of the subsurface from wide‐azimuth seismic data. This method is suitable to build high‐resolution velocity models provided that the targeted area is sampled by both diving waves and reflected waves. However, the conventional formulation of full‐waveform inversion prevents the reconstruction of the small wavenumber components of the velocity model when the subsurface is sampled by reflected waves only. This typically occurs as the depth becomes significant with respect to the length of the receiver array. This study first aims to highlight the limits of the conventional form of full‐waveform inversion when applied to seismic reflection data, through a simple canonical example of seismic imaging and to propose a new inversion workflow that overcomes these limitations. The governing idea is to decompose the subsurface model as a background part, which we seek to update and a singular part that corresponds to some prior knowledge of the reflectivity. Forcing this scale uncoupling in the full‐waveform inversion formalism brings out the transmitted wavepaths that connect the sources and receivers to the reflectors in the sensitivity kernel of the full‐waveform inversion, which is otherwise dominated by the migration impulse responses formed by the correlation of the downgoing direct wavefields coming from the shot and receiver positions. This transmission regime makes full‐waveform inversion amenable to the update of the long‐to‐intermediate wavelengths of the background model from the wide scattering‐angle information. However, we show that this prior knowledge of the reflectivity does not prevent the use of a suitable misfit measurement based on cross‐correlation, to avoid cycle‐skipping issues as well as a suitable inversion domain as the pseudo‐depth domain that allows us to preserve the invariant property of the zero‐offset time. This latter feature is useful to avoid updating the reflectivity information at each non‐linear iteration of the full‐waveform inversion, hence considerably reducing the computational cost of the entire workflow. Prior information of the reflectivity in the full‐waveform inversion formalism, a robust misfit function that prevents cycle‐skipping issues and a suitable inversion domain that preserves the seismic invariant are the three key ingredients that should ensure well‐posedness and computational efficiency of full‐waveform inversion algorithms for seismic reflection data.     

Accelerating the T‐matrix approach to seismic full‐waveform inversion by domain decomposition

A seismic variant of the distorted Born iterative inversion method, which is commonly used in electromagnetic and acoustic (medical) imaging, has been recently developed on the basis of the T‐matrix approach of multiple scattering theory. The distorted Born iterative method is consistent with the Gauss–Newton method, but its implementation is different, and there are potentially significant computational advantages of using the T‐matrix approach in this context. It has been shown that the computational cost associated with the updating of the background medium Green functions after each iteration can be reduced via the use of various linearisation or quasi‐linearisation techniques. However, these techniques for reducing the computational cost may not work well in the presence of strong contrasts. To deal with this, we have now developed a domain decomposition method, which allows one to decompose the seismic velocity model into an arbitrary number of heterogeneous domains that can be treated separately and in parallel. The new domain decomposition method is based on the concept of a scattering‐path matrix, which is well known in solid‐state physics. If the seismic model consists of different domains that are well separated (e.g., different reservoirs within a sedimentary basin), then the scattering‐path matrix formulation can be used to derive approximations that are sufficiently accurate but far more speedy and much less memory demanding because they ignore the interaction between different domains. However, we show here that one can also use the scattering‐path matrix formulation to calculate the overall T‐matrix for a large model exactly without any approximations at a computational cost that is significantly smaller than the cost associated with an exact formal matrix inversion solution. This is because we have derived exact analytical results for the special case of two interacting domains and combined them with Strassen's formulas for fast recursive matrix inversion. To illustrate the fact that we have accelerated the T‐matrix approach to full‐waveform inversion by domain decomposition, we perform a series of numerical experiments based on synthetic data associated with a complex salt model and a simpler two‐dimensional model that can be naturally decomposed into separate upper and lower domains. If the domain decomposition method is combined with an additional layer of multi‐scale regularisation (based on spatial smoothing of the sensitivity matrix and the data residual vector along the receiver line) beyond standard sequential frequency inversion, then one apparently can also obtain stable inversion results in the absence of ultra‐low frequencies and reduced computation times.     

Offset‐variable density improves acoustic full‐waveform inversion: a shallow marine case study

We have previously applied three‐dimensional acoustic, anisotropic, full‐waveform inversion to a shallow‐water, wide‐angle, ocean‐bottom‐cable dataset to obtain a high‐resolution velocity model. This velocity model produced an improved match between synthetic and field data, better flattening of common‐image gathers, a closer fit to well logs, and an improvement in the pre‐stack depth‐migrated image. Nevertheless, close examination reveals that there is a systematic mismatch between the observed and predicted data from this full‐waveform inversion model, with the predicted data being consistently delayed in time. We demonstrate that this mismatch cannot be produced by systematic errors in the starting model, by errors in the assumed source wavelet, by incomplete convergence, or by the use of an insufficiently fine finite‐difference mesh. Throughout these tests, the mismatch is remarkably robust with the significant exception that we do not see an analogous mismatch when inverting synthetic acoustic data. We suspect therefore that the mismatch arises because of inadequacies in the physics that are used during inversion. For ocean‐bottom‐cable data in shallow water at low frequency, apparent observed arrival times, in wide‐angle turning‐ray data, result from the characteristics of the detailed interference pattern between primary refractions, surface ghosts, and a large suite of wide‐angle multiple reflected and/or multiple refracted arrivals. In these circumstances, the dynamics of individual arrivals can strongly influence the apparent arrival times of the resultant compound waveforms. In acoustic full‐waveform inversion, we do not normally know the density of the seabed, and we do not properly account for finite shear velocity, finite attenuation, and fine‐scale anisotropy variation, all of which can influence the relative amplitudes of different interfering arrivals, which in their turn influence the apparent kinematics. Here, we demonstrate that the introduction of a non‐physical offset‐variable water density during acoustic full‐waveform inversion of this ocean‐bottom‐cable field dataset can compensate efficiently and heuristically for these inaccuracies. This approach improves the travel‐time match and consequently increases both the accuracy and resolution of the final velocity model that is obtained using purely acoustic full‐waveform inversion at minimal additional cost.     

1D elastic full‐waveform inversion and uncertainty estimation by means of a hybrid genetic algorithm–Gibbs sampler approach

Stochastic optimization methods, such as genetic algorithms, search for the global minimum of the misfit function within a given parameter range and do not require any calculation of the gradients of the misfit surfaces. More importantly, these methods collect a series of models and associated likelihoods that can be used to estimate the posterior probability distribution. However, because genetic algorithms are not a Markov chain Monte Carlo method, the direct use of the genetic‐algorithm‐sampled models and their associated likelihoods produce a biased estimation of the posterior probability distribution. In contrast, Markov chain Monte Carlo methods, such as the Metropolis–Hastings and Gibbs sampler, provide accurate posterior probability distributions but at considerable computational cost. In this paper, we use a hybrid method that combines the speed of a genetic algorithm to find an optimal solution and the accuracy of a Gibbs sampler to obtain a reliable estimation of the posterior probability distributions. First, we test this method on an analytical function and show that the genetic algorithm method cannot recover the true probability distributions and that it tends to underestimate the true uncertainties. Conversely, combining the genetic algorithm optimization with a Gibbs sampler step enables us to recover the true posterior probability distributions. Then, we demonstrate the applicability of this hybrid method by performing one‐dimensional elastic full‐waveform inversions on synthetic and field data. We also discuss how an appropriate genetic algorithm implementation is essential to attenuate the “genetic drift” effect and to maximize the exploration of the model space. In fact, a wide and efficient exploration of the model space is important not only to avoid entrapment in local minima during the genetic algorithm optimization but also to ensure a reliable estimation of the posterior probability distributions in the subsequent Gibbs sampler step.     

Anisotropic elastic full‐waveform inversion of walkaway vertical seismic profiling data from the Arabian Gulf

Borehole seismic addresses the need for high‐resolution images and elastic parameters of the subsurface. Full‐waveform inversion of vertical seismic profile data is a promising technology with the potential to recover quantitative information about elastic properties of the medium. Full‐waveform inversion has the capability to process the entire wavefield and to address the wave propagation effects contained in the borehole data—multi‐component measurements; anisotropic effects; compressional and shear waves; and transmitted, converted, and reflected waves and multiples. Full‐waveform inversion, therefore, has the potential to provide a more accurate result compared with conventional processing methods. We present a feasibility study with results of the application of high‐frequency (up to 60 Hz) anisotropic elastic full‐waveform inversion to a walkaway vertical seismic profile data from the Arabian Gulf. Full‐waveform inversion has reproduced the majority of the wave events and recovered a geologically plausible layered model with physically meaningful values of the medium.     

Acoustic full waveform inversion of synthetic land and marine data in the Laplace domain

Elastic waves, such as Rayleigh and mode‐converted waves, together with amplitude versus offset variations, serve as noise in full waveform inversion using the acoustic approximation. Heavy preprocessing must be applied to remove elastic effects to invert land or marine data using the acoustic inversion method in the time or frequency domains. Full waveform inversion using the elastic wave equation should be one alternative; however, multi‐parameter inversion is expensive and sensitive to the starting velocity model. We implement full acoustic waveform inversion of synthetic land and marine data in the Laplace domain with minimum preprocessing (i.e., muting) to remove elastic effects. The damping in the Laplace transform can be thought of as an automatic time windowing. Numerical examples show that Laplace‐domain acoustic inversion can yield correct smooth velocity models even with the noise originating from elastic waves. This offers the opportunity to develop an accurate smooth starting model for subsequent inversion in the frequency domain.     

Determining the focal mechanisms and depths of relatively small earthquakes using a few stations by full‐waveform modelling

Determining the focal mechanism of earthquakes helps us to better define faults and understand the stress regime. This technique can be helpful in the oil and gas industry where it can be applied to microseismic events. The objective of this paper is to find double couple focal mechanisms, excluding scalar seismic moments, and the depths of small earthquakes using data from relatively few local stations. This objective is met by generating three‐component synthetic seismograms to match the observed normalized velocity seismograms. We first calculate Green's functions given an initial estimate of the earthquake's hypocentre, the locations of the seismic recording stations and a 1D velocity model of the region for a series of depths. Then, we calculate the moment tensor for different combinations of strikes, dips and rakes for each depth. These moment tensors are combined with the Green's functions and then convolved with a source time function to produce synthetic seismograms. We use a grid search to find the synthetic seismogram with the largest objective function that best fits all three components of the observed velocity seismogram. These parameters define the focal mechanism solution of an earthquake. We tested the method using three earthquakes in Southern California with moment magnitudes of 5.0, 5.1 and 4.4 using the frequency range 0.1–2.0 Hz. The source mechanisms of the events were determined independently using data from a multitude of stations. Our results obtained, from as few as three stations, generally match those obtained by the Southern California Earthquake Data Center. The main advantage of this method is that we use relatively high‐frequency full‐waveforms, including those from short‐period instruments, which makes it possible to find the focal mechanism and depth of earthquakes using as few as three stations when the velocity structure is known.     

基于地震数据子集的波形反演思路、方法与应用

地震数据与地下介质物性参数之间的复杂关系,决定了地震全波形反演在理论方法上面临着强烈的非线性难题.地下不同物性参数的不同分量在地震数据上具有不同的表现,勘探的不同阶段对地下介质模型的精度也具有不同的要求,这就决定了在地震全波形反演过程中不必时刻追求地震数据全部信息的匹配,部分信息的匹配就有可能解决现阶段的某些问题,还可以一定程度上规避匹配全部地震信息所遇到的强烈非线性难题.基于这样的考虑,我们提出了利用地震数据子集进行波形反演的思路,给出了统一的反演方法,并通过基于包络数据子集以及反射波数据子集的波形反演的理论模型与实际资料反演试验,证明了所提出的波形反演思路和方法的正确性.     

Frequency‐domain waveform modelling and inversion for coupled media using a symmetric impedance matrix

To simulate the seismic signals that are obtained in a marine environment, a coupled system of both acoustic and elastic wave equations is solved. The acoustic wave equation for the fluid region simulates the pressure field while minimizing the number of degrees of freedom of the impedance matrix, and the elastic wave equation for the solid region simulates several elastic events, such as shear waves and surface waves. Moreover, by combining this coupled approach with the waveform inversion technique, the elastic properties of the earth can be inverted using the pressure data obtained from the acoustic region. However, in contrast to the pure acoustic and elastic cases, the complex impedance matrix for the coupled media does not have a symmetric form because of the boundary (continuity) condition at the interface between the acoustic and elastic elements. In this study, we propose a manipulation scheme that makes the complex impedance matrix for acoustic–elastic coupled media to take a symmetric form. Using the proposed symmetric matrix, forward and backward wavefields are identical to those generated by the conventional approach; thus, we do not lose any accuracy in the waveform inversion results. However, to solve the modified symmetric matrix, LDLT factorization is used instead of LU factorization for a matrix of the same size; this method can mitigate issues related to severe memory insufficiency and long computation times, particularly for large‐scale problems.     

Optimized experimental design for seismic full waveform inversion: A computationally efficient method including a flexible implementation of acquisition costs

Optimized experimental design aims at reducing the cost of a seismic survey by identifying the optimal locations and amounts of sources and receivers. While the acquisition design in the context of seismic imaging applies criteria like fold, offset and spatial sampling, different attributes such as the sensitivity kernels are more relevant for seismic full waveform inversion. An ideal measure to quantify the goodness of an acquisition design relies on the eigenvalue spectrum of the approximate Hessian matrix, but this technique is computationally too expensive for practical use. A more affordable goodness measure has been proposed in the past, but we demonstrate that this measure is inappropriate for target‐oriented optimized experimental design. To address those issues, we derived a sequential receiver‐based procedure using a goodness measure based on the determinant of the approximate Hessian matrix. We show with numerical tests that it efficiently provides an optimized design for target‐oriented as well as for extensive full waveform inversion. This design allows a better reconstruction of the subsurface than an evenly spaced acquisition geometry. Furthermore, the optimization algorithm itself can easily be parallelized, therefore making it attractive for applications to large‐scale three‐dimensional surveys. In addition, our algorithm is able to incorporate variable costs, representing any kind of acquisition‐related costs, for every individual source location. The combined optimization with respect to the information content of sources and to the true cost will allow a more comprehensive and realistic survey planning and has a high potential for further applications.     

基于多网格的频率域全波形反演(英文)

频率域全波形反演虽然克服了时间方向上的局部极小值问题,但是地下介质的复杂性使其在空间域仍然存在局部极小值缺陷。在优化梯度法基础上,本文采用预条件双共轭梯度稳定算法和多重网格方法计算反演中的波场传播和目标函数的梯度,在保证计算速度的同时,减小计算机内存的消耗。频率域波形反演和多重网格的多尺度性质有效改善问题极小值缺陷,加快反演的收敛速度。以局部非均匀的三孔模型和Marmousi模型的数值模拟结果验证了该算法的有效性。     

Full waveform inversion of SH‐ and Love‐wave data in near‐surface prospecting

We develop a two‐dimensional full waveform inversion approach for the simultaneous determination of S‐wave velocity and density models from SH ‐ and Love‐wave data. We illustrate the advantages of the SH/Love full waveform inversion with a simple synthetic example and demonstrate the method's applicability to a near‐surface dataset, recorded in the village ?achtice in Northwestern Slovakia. Goal of the survey was to map remains of historical building foundations in a highly heterogeneous subsurface. The seismic survey comprises two parallel SH‐profiles with maximum offsets of 24 m and covers a frequency range from 5 Hz to 80 Hz with high signal‐to‐noise ratio well suited for full waveform inversion. Using the Wiechert–Herglotz method, we determined a one‐dimensional gradient velocity model as a starting model for full waveform inversion. The two‐dimensional waveform inversion approach uses the global correlation norm as objective function in combination with a sequential inversion of low‐pass filtered field data. This mitigates the non‐linearity of the multi‐parameter inverse problem. Test computations show that the influence of visco‐elastic effects on the waveform inversion result is rather small. Further tests using a mono‐parameter shear modulus inversion reveal that the inversion of the density model has no significant impact on the final data fit. The final full waveform inversion S‐wave velocity and density models show a prominent low‐velocity weathering layer. Below this layer, the subsurface is highly heterogeneous. Minimum anomaly sizes correspond to approximately half of the dominant Love‐wavelength. The results demonstrate the ability of two‐dimensional SH waveform inversion to image shallow small‐scale soil structure. However, they do not show any evidence of foundation walls.     

Analysis of different parameterisations of waveform inversion of compressional body waves in an elastic transverse isotropic Earth with a vertical axis of symmetry

In a multi‐parameter waveform inversion, the choice of the parameterisation influences the results and their interpretations because leakages and the tradeoff between parameters can cause artefacts. We review the parameterisation selection when the inversion focuses on the recovery of the intermediate‐to‐long wavenumbers of the compressional velocities from the compressional body (P) waves. Assuming a transverse isotropic medium with a vertical axis of symmetry and weak anisotropy, analytical formulas for the radiation patterns are developed to quantify the tradeoff between the shear velocity and the anisotropic parameters and the effects of setting to zero the shear velocity in the acoustic approach. Because, in an anisotropic medium, the radiation patterns depend on the angle of the incident wave with respect to the vertical axis, two particular patterns are discussed: a transmission pattern when the ingoing and outgoing slowness vectors are parallel and a reflection pattern when the ingoing and outgoing slowness vectors satisfy Snell's law. When the inversion aims at recovering the long‐to‐intermediate wavenumbers of the compressional velocities from the P‐waves, we propose to base the parameterisation choice on the transmission patterns. Since the P‐wave events in surface seismic data do not constrain the background (smooth) vertical velocity due to the velocity/depth ambiguity, the preferred parameterisation contains a parameter that has a transmission pattern concentrated along the vertical axis. This parameter can be fixed during the inversion which reduces the size of the model space. The review of several parameterisations shows that the vertical velocity, the Thomsen parameter δ, or the Thomsen parameter ε have a transmission pattern along the vertical axis depending on the parameterisation choice. The review of the reflection patterns of those selected parameterisations should be done in the elastic context. Indeed, when reflection data are also inverted, there are potential leakages of the shear parameter at intermediate angles when we carry out acoustic inversion.     

基于精确震源函数的解调包络多尺度全波形反演

本文提出解调包络方法来重构地震记录中缺失的低频信号,同时该方法能够降低全波形反演的非线性程度;提出伴随状态震源函数反演方法来得到精确的震源函数,并推导了梯度计算公式;解调包络方法结合低通滤波技术,实现了从低频到高频的多尺度反演策略,有效缓解了全波形反演的周波跳跃问题.数值算例证明了解调包络、伴随状态震源函数反演方法和低通滤波多尺度反演策略的可行性及优越性.震源函数反演精度测试结果表明：即使观测记录在缺失低频信息的情况下,也能反演得到精确的震源函数.缺失低频测试和抗噪能力测试结果表明：即使地震数据中缺失9Hz以下的低频信号或者信噪比极低的情况下,利用反演得到的精确震源函数进行解调包络多尺度全波形反演,同样可以得到高精度的全波形反演结果.与Hilbert包络全波形反演对比结果表明：解调包络在重构低频和降低伴随震源主频方面具有一定优势.

     

Improved hybrid iterative optimization method for seismic full waveform inversion

In full waveform inversion (FWI), Hessian information of the misfit function is of vital importance for accelerating the convergence of the inversion; however, it usually is not feasible to directly calculate the Hessian matrix and its inverse. Although the limited memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) or Hessian-free inexact Newton (HFN) methods are able to use approximate Hessian information, the information they collect is limited. The two methods can be interlaced because they are able to provide Hessian information for each other; however, the performance of the hybrid iterative method is dependent on the effective switch between the two methods. We have designed a new scheme to realize the dynamic switch between the two methods based on the decrease ratio (DR) of the misfit function (objective function), and we propose a modified hybrid iterative optimization method. In the new scheme, we compare the DR of the two methods for a given computational cost, and choose the method with a faster DR. Using these steps, the modified method always implements the most efficient method. The results of Marmousi and over thrust model testings indicate that the convergence with our modified method is significantly faster than that in the L-BFGS method with no loss of inversion quality. Moreover, our modified outperforms the enriched method by a little speedup of the convergence. It also exhibits better efficiency than the HFN method.     

地震密度反演及地层孔隙度估计

地层密度直接与孔隙度、孔隙流体类型、饱和度和骨架矿物成分有关.本文通过理论分析和计算,讨论了油气藏储层物性参数变化引起的密度变化及密度变化对地震波速度、阻抗和振幅的影响,提出了基于完全纵波方程的全波形地震密度反演和孔隙度估计方法,克服了常规地震密度反演对地震数据更多处理引起的信号畸变,提高了地震密度反演和地层孔隙度估计的精度.该方法采用波场导数的时间积分和多炮求和,对地震数据中的噪声具有比较强的压制作用.理论模型研究表明该方法是可行的.通过对我国西部某气田实际数据处理、分析和反演,获得了地层密度和孔隙度,结果与测井基本吻合,证明了预测结果的准确性和方法的有效性,从而为后续的有效储层预测和储量计算提供了可靠的数据.     

Correlation‐based reflection waveform inversion by one‐way wave equations

Reflection full waveform inversion can update subsurface velocity structure of the deeper part, but tends to get stuck in the local minima associated with the waveform misfit function. These local minima cause cycle skipping if the initial background velocity model is far from the true model. Since conventional reflection full waveform inversion using two‐way wave equation in time domain is computationally expensive and consumes a large amount of memory, we implement a correlation‐based reflection waveform inversion using one‐way wave equations to retrieve the background velocity. In this method, one‐way wave equations are used for the seismic wave forward modelling, migration/de‐migration and the gradient computation of objective function in frequency domain. Compared with the method using two‐way wave equation, the proposed method benefits from the lower computational cost of one‐way wave equations without significant accuracy reduction in the cases without steep dips. It also largely reduces the memory requirement by an order of magnitude than implementation using two‐way wave equation both for two‐ and three‐dimensional situations. Through numerical analysis, we also find that one‐way wave equations can better construct the low wavenumber reflection wavepath without producing high‐amplitude short‐wavelength components near the image points in the reflection full waveform inversion gradient. Synthetic test and real data application show that the proposed method efficiently updates the background velocity model.     

基于声学全波形反演的油气藏地震成像方法

岩性油气藏在我国天然气勘探开发中占有非常重要的位置,其分布区域的成像是合理布设井位,提高钻井成功率的关键之一.本文首先基于地下介质的声学近似和波场回传理论,利用频率域单程声波方程延拓计算地震波场,进行全波形反演,获得地层密度和体积模量的定量成像,并依据油气藏物性特征和流体饱和多孔介质岩石物理模型,简要讨论了孔隙度和饱和度与密度及体积模量的关系,明确了地震油气藏成像新概念.在此基础上,定义了基于流体体积模量和孔隙度的成像函数,进行油气藏成像.理论模型计算表明该方法是可行的.通过对西部地区某气田二维地震数据处理,实现了致密砂岩气藏成像,钻井结果证实了气藏区域成像位置的准确性和方法的有效性.     

Sensitivity analysis and application of time‐lapse full‐waveform inversion: synthetic testing and field data example from the North Sea,Norway

Time‐lapse refraction can provide complementary seismic solutions for monitoring subtle subsurface changes that are challenging for conventional P‐wave reflection methods. The utilization of refraction time lapse has lagged behind in the past partly due to the lack of robust techniques that allow extracting easy‐to‐interpret reservoir information. However, with the recent emergence of the full‐waveform inversion technique as a more standard tool, we find it to be a promising platform for incorporating head waves and diving waves into the time‐lapse framework. Here we investigate the sensitivity of 2D acoustic, time‐domain, full‐waveform inversion for monitoring a shallow, weak velocity change (?30 m/s, or ?1.6%). The sensitivity tests are designed to address questions related to the feasibility and accuracy of full‐waveform inversion results for monitoring the field case of an underground gas blowout that occurred in the North Sea. The blowout caused the gas to migrate both vertically and horizontally into several shallow sand layers. Some of the shallow gas anomalies were not clearly detected by conventional 4D reflection methods (i.e., time shifts and amplitude difference) due to low 4D signal‐to‐noise ratio and weak velocity change. On the other hand, full‐waveform inversion sensitivity analysis showed that it is possible to detect the weak velocity change with the non‐optimal seismic input. Detectability was qualitative with variable degrees of accuracy depending on different inversion parameters. We inverted, the real 2D seismic data from the North Sea with a greater emphasis on refracted and diving waves’ energy (i.e., most of the reflected energy was removed for the shallow zone of interest after removing traces with offset less than 300 m). The full‐waveform inversion results provided more superior detectability compared with the conventional 4D stacked reflection difference method for a weak shallow gas anomaly (320 m deep).     

Determination of porosity and saturation using seismic waveform inversion

Characterization of a reservoir model requires determination of its petrophysical parameters, such as porosity and saturation. We propose a new method to determine these parameters directly from seismic data. The method consists of the computation and inversion of seismic waveforms. A high frequency method is presented to model wave propagation through an attenuative and dispersive poroelastic medium. The high frequency approximation makes it possible to efficiently compute sensitivity functions. This enables the inversion of seismic waveforms for porosity and saturation. The waveform inversion algorithm is applied to two laboratory crosswell datasets of a water saturated sand. The starting models were obtained using travel time tomography. The first dataset is inverted for porosity. The misfit reduction for this dataset is approximately 50%. The second dataset was obtained after injection of a nonaqueous-phase liquid (NAPL), possibly with some air, which made the medium more heterogeneous. This dataset was inverted for NAPL and air saturation using the porosity model obtained from the first inversion. The misfit reduction of the second experiment was 70%. Regions of high NAPL and high air saturation were found at the same location. These areas correlate well with the position of one of the injection points as well as regions of higher NAPL concentrations found after excavation of the sand. It is therefore possible to directly invert waveforms for pore fluid saturation by taking into account the attenuation and dispersion caused by the poroelasticity.