Marine vibrators: the new phase of seismic exploration

Marine seismic vibrators are generally considered to be less intrusive than airguns from an environmental perspective. This is because they emit their energy spread out in time, rather than in a single, high-intensity pulse. There are also significant geophysical benefits associated with marine vibrators, and they stem from the ability to specify in detail the output acoustic waveform. The phase can be specified independently at each frequency. Such detailed control cannot be achieved with conventional airgun sources, where the phase can only be modified using simple overall time delays. The vibrator phase can be employed in several different ways: it can be applied to the overall source phase in a sequence so that it varies from one source point to the next; it can be applied to the individual vibrators within the source array so the source directivity is changed; it can be applied to the overall source phase of each source in a simultaneous source acquisition. Carefully designed phase sequences can attenuate the residual source noise, and this in turn allows extra source points to be interleaved between the conventional ones. For these extra source points, the relative phase of the vibrators within the array can be chosen to create a transverse gradient source, which illuminates the earth predominantly in directions out of the plane of the sail line without left/right ambiguity. If seismic vibrator data are acquired using interleaved conventional and transverse gradient sweeps, more information is collected per kilometre of vessel travel than is the case in conventional acquisition. This richer data acquisition leads to the possibility of acquiring all the necessary seismic data in a shorter time. Three-dimensional reconstruction techniques are used to recover the same image quality that would have been obtained using the conventional, more time-consuming acquisition. For a marine vibrator to be suitable for these techniques it must, in general terms, have ‘high fidelity’. The precise device specifications are defined through realistic end-to-end simulations of the physical systems and the processing. The specifications are somewhat more onerous than for a conventional vibrator, but they are achievable. A prototype vibrator that satisfies these requirements has been built. In a simulated case study of a three-dimensional deep-water ocean bottom node survey, the seismic data could have been acquired using marine vibrators in one third of the time that it would have taken using airguns.     

Interaction of helically wound fibre‐optic cables with plane seismic waves

Distributed acoustic sensing is a novel technology for seismic acquisition. In this technology, strain changes induced by seismic waves impinging on an optical fibre are monitored. Due to the fact that glass is relatively rigid, straight glass fibres are not sensitive to broadside waves. We suggest using distributed acoustic sensing systems with fibres helically wound around cables. One increases the fibre sensitivity to broadside waves by decreasing the fibre wrapping angle (the angle between the fibre axis and the plane normal to the cable axis). The optimal wrapping angle is chosen to minimize the impact of Rayleigh waves on the signal measured. This angle depends on the cable Poisson ratio, and it is approximately equal to 30° for cables composed of plastic. For reliable detection of seismic waves, one needs a good mechanical contact between the cable and the surrounding medium. On the other hand, the sensitivity of distributed acoustic sensing systems to primary waves can be significantly reduced if the cable is placed in a cemented borehole.     

CO2 storage in the high Arctic: efficient modelling of pre‐stack depth‐migrated seismic sections for survey planning

The sequestration of CO₂ in subsurface reservoirs constitutes an immediate counter‐measure to reduce anthropogenic emissions of CO₂, now recognized by international scientific panels to be the single most critical factor driving the observed global climatic warming. To ensure and verify the safe geological containment of CO₂ underground, monitoring of the CO₂ site is critical. In the high Arctic, environmental considerations are paramount and human impact through, for instance, active seismic surveys, has to be minimized. Efficient seismic modelling is a powerful tool to test the detectability and imaging capability prior to acquisition and thus improve the characterization of CO₂ storage sites, taking both geological setting and seismic acquisition set‐up into account. The unique method presented here avoids the costly generation of large synthetic data sets by employing point spread functions to directly generate pre‐stack depth‐migrated seismic images. We test both a local‐target approach using an analytical filter assuming an average velocity and a full‐field approach accounting for the spatial variability of point spread functions. We assume a hypothetical CO₂ plume emplaced in a sloping aquifer inspired by the conditions found at the University of Svalbard CO₂ lab close to Longyearbyen, Svalbard, Norway, constituting an unconventional reservoir–cap rock system. Using the local‐target approach, we find that even the low‐to‐moderate values of porosity (5%–18%) measured in the reservoir should be sufficient to induce significant change in seismic response when CO₂ is injected. The sensitivity of the seismic response to changes in CO₂ saturation, however, is limited once a relatively low saturation threshold of 5% is exceeded. Depending on the illumination angle provided by the seismic survey, the quality of the images of five hypothetical CO₂ plumes of varying volume differs depending on the steepness of their flanks. When comparing the resolution of two orthogonal 2D surveys to a 3D survey, we discover that the images of the 2D surveys contain significant artefacts, the CO₂‐brine contact is misplaced and an additional reflector is introduced due to the projection of the point spread function of the unresolvable plane onto the imaging plane. All of these could easily lead to a misinterpretation of the behaviour of the injected CO₂. Our workflow allows for testing the influence of geological heterogeneities in the target aquifer (igneous intrusions, faults, pervasive fracture networks) by utilizing increasingly complex and more realistic geological models as input as more information on the subsurface becomes available.     

A study of the geophysical response of distributed fibre optic acoustic sensors through laboratory‐scale experiments

In the past few years, distributed acoustic sensing has gained great interest in geophysics. This acquisition technology offers immense improvement in terms of efficiency when compared with current geophysical acquisition methods. However, the fundamentals of the measurement are still not fully understood because direct comparisons of fibre data with conventional geophysical sensors are difficult during field tests. We present downscaled laboratory experiments that enabled us to characterise the relationship between the signals recorded by conventional seismic point receivers and by distributed fibre optic sensors. Interrogation of the distributed optical fibre sensor was performed with a Michelson interferometer because this system is suited to compact test configurations, and it requires only a very simple data processing workflow for extracting the signal outputs. We show acoustic data that were recorded simultaneously by both the fibre optical interferometer and conventional three‐component accelerometers, thus enabling the comparison of sensor performance. We present results focused on the directionality of fibre measurements, on the amplitude variation with angle of incidence, and on the transfer function that allows accelerometer signals to be transformed into optical fibre signals. We conclude that the optical fibre response matches with the array of the displacement differences of the inline accelerometers deployed along the fibre length. Moreover, we also analysed the influence of various types of coupling and fibre cable coating on the signal responses, emphasising the importance of these parameters for field seismic acquisitions when using the distributed fibre optic technology.     

Field testing of modular borehole monitoring with simultaneous distributed acoustic sensing and geophone vertical seismic profiles at Citronelle,Alabama

A modular borehole monitoring concept has been implemented to provide a suite of well‐based monitoring tools that can be deployed cost effectively in a flexible and robust package. The initial modular borehole monitoring system was deployed as part of a CO₂ injection test operated by the Southeast Regional Carbon Sequestration Partnership near Citronelle, Alabama. The Citronelle modular monitoring system transmits electrical power and signals, fibre‐optic light pulses, and fluids between the surface and a reservoir. Additionally, a separate multi‐conductor tubing‐encapsulated line was used for borehole geophones, including a specialized clamp for casing clamping with tubing deployment. The deployment of geophones and fibre‐optic cables allowed comparison testing of distributed acoustic sensing. We designed a large source effort (>64 sweeps per source point) to test fibre‐optic vertical seismic profile and acquired data in 2013. The native measurement in the specific distributed acoustic sensing unit used (an iDAS from Silixa Ltd) is described as a localized strain rate. Following a processing flow of adaptive noise reduction and rebalancing the signal to dimensionless strain, improvement from repeated stacking of the source was observed. Conversion of the rebalanced strain signal to equivalent velocity units, via a scaling by local apparent velocity, allows quantitative comparison of distributed acoustic sensing and geophone data in units of velocity. We see a very good match of uncorrelated time series in both amplitude and phase, demonstrating that velocity‐converted distributed acoustic sensing data can be analyzed equivalent to vertical geophones. We show that distributed acoustic sensing data, when averaged over an interval comparable to typical geophone spacing, can obtain signal‐to‐noise ratios of 18 dB to 24 dB below clamped geophones, a result that is variable with noise spectral amplitude because the noise characteristics are not identical. With vertical seismic profile processing, we demonstrate the effectiveness of downgoing deconvolution from the large spatial sampling of distributed acoustic sensing data, along with improved upgoing reflection quality. We conclude that the extra source effort currently needed for tubing‐deployed distributed acoustic sensing vertical seismic profile, as part of a modular monitoring system, is well compensated by the extra spatial sampling and lower deployment cost as compared with conventional borehole geophones.     

Target-level waveform inversion: a prospective application of the convolution-type representation for the acoustic wavefield

Nowadays, full-waveform inversion, based on fitting the measured surface data with modelled data, has become the preferred approach to recover detailed physical parameters from the subsurface. However, its application is computationally expensive for large inversion domains. Furthermore, when the subsurface has a complex geological setting, the inversion process requires an appropriate pre-conditioning scheme to retrieve the medium parameters for the desired target area in a reliable manner. One way of dealing with both aspects is by waveform inversion schemes in a target-oriented fashion. Therefore, we propose a prospective application of the convolution-type representation for the acoustic wavefield in the frequency–space domain formulated as a target-oriented waveform inversion method. Our approach aims at matching the observed and modelled upgoing wavefields at a target depth level in the subsurface, where the seismic wavefields, generated by sources distributed above this level, are available. The forward modelling is performed by combining the convolution-type representation for the acoustic wavefield with solving the two-way acoustic wave-equation in the frequency–space domain for the target area. We evaluate the effectiveness of our inversion method by comparing it with the full-domain full-waveform inversion process through some numerical examples using synthetic data from a horizontal well acquisition geometry, where the sources are located at the surface and the receivers are located along a horizontal well at the target level. Our proposed inversion method requires less computational effort and, for this particular acquisition, it has proven to provide more accurate estimates of the target zone below a complex overburden compared to both full-domain full-waveform inversion process and local full-waveform inversion after applying interferometry by multidimensional deconvolution to get local-impulse responses.     

A genetic algorithm for filter design to enhance features in seismic images

We present a novel method to enhance seismic data for manual and automatic interpretation. We use a genetic algorithm to optimize a kernel that, when convolved with the seismic image, appears to enhance the internal characteristics of salt bodies and the sub‐salt stratigraphy. The performance of the genetic algorithm was validated by the use of test images prior to its application on the seismic data. We present the evolution of the resulting kernel and its convolved image. This image was analysed by a seismic interpreter, highlighting possible advantages over the original one. The effects of the kernel were also subject to an automatic interpretation technique based on principal component analysis. Statistical comparison of these results with those from the original image, by means of the Mann‐Whitney U‐test, proved the convolved image to be more appropriate for automatic interpretation.     

Application of the redundant-lifting scheme for ground-roll attenuation in near-surface characterization using full-waveform inversion on P-wave seismic data

Seismic modelling of the shallow subsurface (within the first few metres) is often challenging when the data are dominated by ground-roll and devoid of reflection. We showed that, even when transmission is the only available phase for analysis, fine-scale and interpretable P-wave velocity (V_P) and attenuation (Q_P⁻¹) models can still be prepared using full-waveform inversion, with data being preconditioned for ground-roll. To prove this idea, we suppressed the ground-roll in two different ways before full-waveform inversion modelling: first, through a bottom mute; second, through a novel wavelet transform-based method known as the redundant-lifting scheme. The applicability of full-waveform inversion is tested through imaging two buried targets. These include a pair of utility water pipes with known diameters of 0.8 m and burial depths of 1.5 m, respectively. The second target is the poorly documented backfill, which was the former location of the pipe(s). The data for full-waveform inversion are acquired along a 2D profile using a static array of 24, 40 Hz vertical component geophones and a buried point source. The results show that (a) the redundant-lifting scheme better suppresses the ground roll, which in turn provides better images of the targets in full-waveform inversion; and (b) the V_P and Q_P⁻¹ models from full-waveform inversion of redundant-lifting scheme data could detect the two targets adequately.     

Rock physics model of tight oil siltstone for seismic prediction of brittleness

Tight oil siltstones are rocks with complex structure at pore scale and are characterized by low porosity and low permeability at macroscale. The production of tight oil siltstone reservoirs can be increased by hydraulic fracturing. For optimal fracking results, it is desirable to map the ability to fracture based on seismic data prior to fracturing. Brittleness is currently thought to be a key parameter for evaluating the ability to fracture. To link seismic information to the brittleness distribution, a rock physics model is required. Currently, there exists no commonly accepted rock physics model for tight oil siltstones. Based on the observed correlation between porosity and mineral composition and known microstructure of tight oil siltstone in Daqing oilfield of Songliao basin, we develop a rock physics model by combining the Voigt–Reuss–Hill average, self-consistent approximation and differential effective medium theory. This rock physics model allows us to explore the dependence of the brittleness on porosity, mineral composition, microcrack volume fraction and microcrack aspect ratio. The results show that, as quartz content increases and feldspar content decreases, Young's modulus tends to increase and Poisson ratio decreases. This is taken as a signature of higher brittleness. Using well log data and seismic inversion results, we demonstrate the versatility of the rock physics template for brittleness prediction.     

A comparison of seismic saltbody interpretation via neural networks at sample and pattern levels

Saltbodies are important subsurface structures that have significant implications for hydrocarbon accumulation and sealing in petroleum reservoirs, and accurate saltbody imaging and delineation is now greatly facilitated with the availability of three-dimensional seismic surveying. However, with the growing demand for larger survey coverage and higher imaging resolution, the size of seismic data is increasing dramatically. Correspondingly, manual saltbody interpretation fails to offer an efficient solution, particularly in exploration areas of complicated salt intrusion history. Recently, artificial intelligence is attracting great attention from geoscientists who desire to utilize the popular machine learning technologies for evolving the interpretational tools capable of mimicking an experienced interpreter's intelligence. This study first implements two popular machine learning tools, the multi-layer perceptron and the convolutional neural network, for delineating seismic saltbodies at sample and pattern levels, respectively, then compares their performance through applications to the synthetic SEAM seismic volume, and moreover tentatively investigates what contributes to the better convolutional neural network delineation. Specifically, the multi-layer perceptron scheme is capable of efficiently utilizing an interpreter's knowledge by selecting, pre-conditioning and integrating a set of seismic attributes that best highlight the target saltbodies, whereas the convolutional neural network scheme makes it possible for saltbody delineation directly from seismic amplitude and thus significantly reduces the dependency on attribute selection from interpreters. It is concluded that the better performance from the convolutional neural network scheme results from two factors. First, the convolutional neural network builds the mapping relationship between the seismic signals and the saltbodies using the original seismic amplitude instead of manually selected seismic attributes, so that the negative impact of using less representative attributes is virtually eliminated. Second and more importantly, the convolutional neural network defines, learns and identifies the saltbodies by utilizing local seismic reflection patterns, so that the seismic noises and processing artefacts of distinct patterns are effectively identified and excluded.     

Optimum size and density of surface grid arrays for retrieving accurate shear-tensile fracturing of microearthquakes

Surface arrays became an important tool for monitoring the induced seismicity in hydraulic fracturing experiments and for assessing the impact of fluid injection on the fracturing process of microearthquakes. The layout of sensors plays a key role in this task because it controls the accuracy of event locations and retrieved seismic moment tensors. We simulate various configurations of grid sensor arrays characterized by a different number of sensors, array span, sensor spacing, depth of sources and various shear/tensile source mechanisms of events. The moment tensor inversion is carried out using synthetically calculated P-wave amplitudes with added random noise. A bias in the solutions is evaluated by errors in the double-couple percentage of inverted moment tensors because the double-couple errors inform us about the sensitivity of the network to detect the shear/tensile fracturing mode of induced microearthquakes. The results show that the accuracy of the double-couple percentage is mostly controlled by the offset-to-depth ratio R defined as the ratio of half of the network size to the event depth. The optimum value of R is in the range of 0.75–1.5 irrespective of the type of the focal mechanism. If 121 (11 × 11) sensors are distributed in a regular grid and recorded data are characterized by a 10% random noise, the double-couple error is less than 6%. This error increases, if R is not optimum or if the number of sensors is reduced. However, even sparse arrays with 49 (7 × 7) or 16 (4 × 4) sensors can yield a reasonable accuracy, provided the surface grid arrays are designed to have an optimum size.     

Source wavelet correction for practical Marchenko imaging: A sub-salt field-data example from the Gulf of Mexico

Imaging a target zone below a salt body can be challenging because large velocity contrasts in the overburden between the salt and surrounding sediments generate internal multiples, which interfere with primary reflections from the target level in the imaging process. This can lead to an erroneous interpretation of reflections in the sub-salt area if multiples are misinterpreted as primaries. The Marchenko redatuming method may enable imaging of the sub-salt target area where the effect of the multiply-scattering overburden is removed. This is achieved by creating a redatumed reflection response where virtual sources and receivers are located below the overburden using a macromodel of the velocity field and the surface reflection data. The accuracy of the redatumed data and the associated internal multiple removal, however, depends on the accurate knowledge of the source wavelet of the acquired reflection data. For the first time, we propose a method which can accurately and reliably correct the amplitudes of the reflection response in field data as required by the Marchenko method. Our method operates by iteratively and automatically updating the source function so as to cancel the most artefact energy in the focusing functions, which are also generated by the Marchenko method. We demonstrate the method on a synthetic dataset and successfully apply it to a field dataset acquired in a deep-water salt environment in the Gulf of Mexico. After the successful source wavelet estimation for the field dataset, we create sub-salt target-oriented images with Marchenko redatumed data. Marchenko images using the proposed source wavelet estimation show clear improvements, such as increased continuity of reflectors, compared to surface-based images and to conventional Marchenko images computed without the inverted source wavelet. Our improvements are corroborated by evidence in the literature and our own synthetic results.     

Improved structural interpretation of legacy 3D seismic data from Karee platinum mine (South Africa) through the application of novel seismic attributes

Seismic detection of faults, dykes, potholes and iron-rich ultramafic pegmatitic bodies is of great importance to the platinum mining industry, as these structures affect safety and efficiency. The application of conventional seismic attributes (such as instantaneous amplitude, phase and frequency) in the hard-rock environment is more challenging than in soft-rock settings because the geology is often complex, reflections disrupted and the seismic energy strongly scattered. We have developed new seismic attributes that sharpen seismic reflections, enabling additional structural information to be extracted from hard-rock seismic data. The symmetry attribute is based on the invariance of an object with respect to transformations such as rotation and reflection; it is independent of the trace reflection amplitude, and hence a better indicator of the lateral continuity of thin and weak reflections. The reflection-continuity detector attribute is based on the Hilbert transform; it enhances the visibility of the peaks and troughs of the seismic traces, and hence the continuity of weak reflections. We demonstrate the effectiveness of these new seismic attributes by applying them to a legacy 3D seismic data set from the Bushveld Complex in South Africa. These seismic attributes show good detection of deep-seated thin (∼1.5 m thick) platinum ore bodies and their associated complex geological structures (faults, dykes, potholes and iron-rich ultramafic pegmatites). They provide a fast, cost-effective and efficient interpretation tool that, when coupled with horizon-based seismic attributes, can reveal structures not seen in conventional interpretations.     

New insights from legacy seismic data: reprocessing of legacy 2D seismic data for imaging of iron-oxide mineralization near Sishen Mine,South Africa

Two overlapping legacy seismic profiles, 130 km long end to end, were shot in the 1990s over the Kuruman Hills on the western margin of the Kaapvaal Craton in southern Africa. The 6-s profiles were aimed at investigating the crustal structure of the western Kaapvaal Craton as well as to locate potential continuation of the Witwatersrand gold-bearing horizons beneath the cover rocks, the latter of which was unsuccessful. In this study, the legacy seismic data are reprocessed and used to image the iron-oxide (mainly haematite) mineralization found in the Kuruman Formation of the Griqualand-West Supergroup, which outcrops along the two seismic profiles. The seismic profiles are located close to the Sishen open pit iron mine, where one of the world's largest iron ore concentrations (986 Mt) is mined. The reprocessed and merged seismic data are combined with magnetic, magnetotelluric, borehole and outcrop data to constrain the interpretation, and all indicate the mineralization host rocks to have ∼500 m thickness and 950 m depth. The seismic data further reveal seismic reflections associated with multiple iron ore horizons, which are affected by a first-order scale syncline and numerous near-vertically dipping (∼65–80°) normal and reverse faults of various orientations and throws, thus providing insight into the structurally controlled iron ore mineralization in the area. Seismic tomography and magnetotellurics characterize the sediments to have a velocity ranging between 5000 and 6000 m/s and a resistivity of <10 Ωm. The seismic imaging of the syncline and associated structural disruptions is important for future mining purposes and plans in the area as these structures might have preserved iron-oxide mineralization from erosion. The reprocessed data thus provide information that could be incorporated in potential future underground mine planning in the area, improving the resource evaluation of the iron-oxide deposit. Legacy seismic data are thus shown to hold intrinsic quality and possible untapped potential that can be realized via data reprocessing.     

间断的Galerkin方法在地震波场数值模拟中的应用概述

通过数值求解描述地震波传播的微分方程来模拟波的传播角度而言,在日趋繁荣的地震波场数值模拟方法的探索中,间断的Galerkin方法(Discontinuous Galerkin,DG)是继有限差分法、伪谱法、有限元法以及谱元法之后出现的另一种有效的数值模拟方法,该方法从2005年以来在地震波场模拟中得到了迅速发展,给这一学科注入了新的活力.但是从国内的发展来看,在地震波场数值模拟方面还没有得到足够的关注.尤其是现有的综述性文献都忽略了这一方法得到的有益的成果,本文将重点介绍间断的Galerkin方法的发展及其在地震波场数值模拟的最新应用研究.     

A methodology to estimate the maximum offset and reflection imaging radius for a 2D/3D vertical seismic profiling survey

We have developed a straightforward and ray based methodology to estimate both the maximum offset and reflection imaging radius for multi‐layered velocity models, which can be used for a 2D/3D VSP survey design. Through numerical examples, we demonstrate that the presence of a high‐velocity layer above a target zone significantly reduces the maximum offset and reflection imaging radius. Our numerical examples also show that including in a migration VSP data acquired beyond a recommended maximum offset, radically degrades the quality of the final VSP image. In addition, unlike the conventional straight‐line based approximation that often produces an incorrect large reflection imaging radius, our methodology predicts the VSP imaging radius with more accuracy than does the conventional approximation.     

Quantitative estimation of water storage and residence time in the epikarst with time‐lapse refraction seismic

The hydrodynamic characterization of the epikarst, the shallow part of the unsaturated zone in karstic systems, has always been challenging for geophysical methods. This work investigates the feasibility of coupling time‐lapse refraction seismic data with petrophysical and hydrologic models for the quantitative determination of water storage and residence time at shallow depth in carbonate rocks. The Biot–Gassmann fluid substitution model describing the seismic velocity variations with water saturation at low frequencies needs to be modified for this lithology. I propose to include a saturation‐dependent rock‐frame weakening to take into account water–rock interactions. A Bayesian inversion workflow is presented to estimate the water content from seismic velocities measured at variable saturations. The procedure is tested first with already published laboratory measurements on core samples, and the results show that it is possible to estimate the water content and its uncertainty. The validated procedure is then applied to a time‐lapse seismic study to locate and quantify seasonal water storage at shallow depth along a seismic profile. The residence time of the water in the shallow layers is estimated by coupling the time‐lapse seismic measurements with rainfall chronicles, simple flow equations, and the petrophysical model. The daily water input computed from the chronicles is used to constraint the inversion of seismic velocities for the daily saturation state and the hydrodynamic parameters of the flow model. The workflow is applied to a real monitoring case, and the results show that the average residence time of the water in the epikarst is generally around three months, but it is only 18 days near an infiltration pathway. During the winter season, the residence times are three times shorter in response to the increase in the effective rainfall.     

射线法模拟分析井间地震观测的波场特征

按照井间地震的观测系统,用改进的突变点加插值射线追踪方法,追踪每炮每道的射线路径,计算几种主要类型的波沿射线路径的波至时间和射线振幅,制作井间地震多炮多道水平分量和垂直分量的合成记录.并将合成记录选排为井间共炮点道集、共接收点道集、共偏移距道集和共中心深度点道集,系统地分析了不同道集内几种主要类型的地震波的传播特征.对野外观测的实际井间地震记录进行了模拟,从复杂的井间地震记录中,识别出井间地震实际观测到的不同类型的波场,为随后的井间地震资料处理和应用提供了依据.     

Mapping of underground cavities by the passive seismic standing waves method: the case study of Barsukovskaya cave (Novosibirsk region,Russia)

This article presents the results of mapping a karst cave by the passive seismic standing waves method. Barsukovskaya cave is located about 100 km southeast of the city of Novosibirsk (Russia). The total length of the cave's passages and grottoes is estimated at about 200 m, the maximum depth from the earth's surface is about 19 m. The method for studying underground cavities used is based on the effect of the generation of standing waves by microtremor in the space between the earth's surface and the cave roof. The accumulation of amplitude spectra of a large number of microtremor records makes it possible to determine the frequencies of the first few modes of these waves. Areal passive seismic survey on the earth's surface above the cave made it possible to construct a map of the lowest mode frequency distribution over the cave roof. Since no standing waves were observed at other points, this map reflects the cave structure in plan, which confirms the comparison with the cave diagram drawn up earlier by one of the speleologists. A schematic map of the depth of the cave roof was constructed using the longitudinal wave velocity V_p = 3120 m/s determined by the rock samples selected near the entrance to the cave. This map at a qualitative level also agrees with the data of speleologists, which indicate that the cave, on average, gradually becomes deeper from the entrance to its dead‐end branches. The shallower depths in comparison with the data of speleologists are apparently explained by a very low estimate of the velocity determined from a rock sample taken near the entrance to the cave. The reliability of the obtained cave mapping results is confirmed by the numerical simulation results using the finite‐element method.