Wolfram Mathematica in Geophysics
Thanks to blog author Anton Ekimenko for his talk
This note was written in the wake of the Wolfram Russian Technology Conference and contains a synopsis of the report I gave. The event took place in June in the city of St. Petersburg. Given the fact that I work in the quarter from the conference venue, I could not help but attend this event. In 2016 and 2017, I listened to the conference reports, and this year I made a presentation. Firstly, an interesting (it seems to me) topic has appeared that we are developing with Kirill Belov , and secondly, after a long study of the legislation of the Russian Federation regarding sanctions policy, two Wolfram Mathematica licenses have appeared in the enterprise where I work.
Before turning to the topic of my speech, I would like to note the good organization of the event. The conference page uses the image of the Kazan Cathedral. The Cathedral is one of the main attractions of St. Petersburg and it is very clearly visible from the hall in which the conference was held.
')
At the entrance to St. Petersburg State University of Economics, participants were met by student assistants - they did not allow them to get lost. During registration, small souvenirs were given out (a toy - flashing adhesions, a pen, stickers with Wolfram symbols). Lunch and coffee breaks were also included in the conference schedule. About delicious coffee and pies, I already noted on the wall of the group - well done chefs. With this introductory part, I would like to emphasize that the event itself, its format and venue already bring positive emotions.
The report, which was prepared by me and Kirill Belov, is called “Using Wolfram Mathematica to solve problems of applied geophysics. Spectral analysis of seismic data or "where the ancient rivers ran." The content of the report covers two parts: firstly, the use of algorithms available in Wolfram Mathematica for the analysis of geophysical data, and secondly, this is how to place geophysical data in Wolfram Mathematica.
First you need to do a little excursion into geophysics. Geophysics is a science that studies the physical properties of rocks. Well, since the rocks have different properties: electrical, magnetic, elastic, then there are appropriate methods of geophysics: electrical exploration, magnetic exploration, seismic exploration ... In the context of this article, we will only touch on seismic exploration. Seismic exploration is the main method for finding oil and gas. The method is based on the excitation of elastic vibrations and the subsequent recording of the response from rocks composing the study area. Excitation of vibrations is carried out on land (dynamite or non-explosive vibration sources of elastic vibrations) or at sea (air guns). Elastic vibrations propagate through the thickness of the rocks, being refracted and reflected at the boundaries of layers with different properties. The reflected waves return to the surface and are recorded by land-based geophones (usually electrodynamic devices based on the movement of a magnet suspended in a coil) or hydrophones in the sea (based on the piezoelectric effect). By the time of arrival of the waves, one can judge the depths of the geological layers.
Seismic vessel tows equipment:
Air gun excites elastic vibrations:
Waves pass through the rock mass and are recorded by hydrophones:
Geophysical exploration research vessel “Ivan Gubkin” on the pier near the Blagoveshchensky bridge in St. Petersburg:
Rocks have different physical properties. For seismic exploration, elastic properties — the propagation velocity of elastic vibrations and density — are of primary importance. If two layers have the same or similar properties, the wave will not “notice” the boundary between them. If the wave velocities in the layers differ, then reflection will occur at the boundary of the layers. The greater the difference in properties, the more intense the reflection. Its intensity will be determined by the reflection coefficient (rc):
where ρ is the rock density, ν is the wave velocity, 1 and 2 denote the upper and lower layers.
One of the simplest and most commonly used seismic signal models is the convolution model, when a registered seismic trace is presented as a result of convolution of a sequence of reflection coefficients with a probe pulse:
where s (t) is the seismic track, i.e. all that a hydrophone or a geophone recorded over a fixed recording time, w (t) is the signal that the airgun generates, n (t) is random noise.
Let us calculate, for example, a synthetic seismic survey. We will use the Ricker pulse widely used in seismic exploration as the initial signal.
Source seismic pulse
We set two boundaries at depths of 300 ms and 600 ms, and the reflection coefficients will be random numbers.
The sequence of reflection coefficients:
Calculate and display the seismic track. Since the reflection coefficients have different signs, we obtain two alternating reflections on the seismic path.
Simulated track:
For this example, a reservation must be made - in reality, the depth of the layers is determined of course in meters, and the calculation of the seismic path occurs for the time domain. It would be more correct to set the depths in meters and calculate the arrival times knowing the velocity in the strata. In this case, I immediately set the layers on the time axis.
If we talk about field studies, as a result of such observations, a huge number of such time series (seismic traces) are recorded. For example, when exploring a site 25 km long and 15 km wide, where as a result of work each track characterizes a cell measuring 25x25 meters (this cell is called a bin), the final data array will contain 600,000 tracks. With a time sampling step of 1 ms, a recording time of 5 seconds, the final data file will be more than 11 GB, and the amount of the original “raw” material can be hundreds of gigabytes.
How to work with them in Wolfram Mathematica ?
The development of the package began with a question on the VK wall of the Russian-speaking support group. Thanks to community responses, a solution was found very quickly. And as a result, it grew into a serious development. The corresponding post on the wall of Wolfram Comunity was even marked by moderators. Currently, the package supports working with the following types of data that are actively used in the geological industry:
To install the package, you must follow the instructions on the download page of the assembled package, i.e. execute the following code on any Mathematica notebook :
After which the package will be installed in the default folder, the path to which can be obtained as follows:
For example, we demonstrate the main features of the package. The call is traditionally made for packages in Wolfram Language:
The package is being developed using Wolfram Workbench . This allows you to accompany the main functionality of the package with documentation, which in presentation format does not differ from the documentation of Wolfram Mathematica itself and provide the package with test files for the first acquaintance.
Such a file, in particular, is the Marmousi.segy file, a synthetic model of a geological section developed by the French Petroleum Institute. Using this model, developers test their own algorithms for modeling the wave field, data processing, inversion of seismic traces, etc. The Marmousi model itself is stored in the repository, from where the package itself was downloaded. In order to get the file, execute the following code:
The result of the import is a SEGYData object:
The SEGY format involves the storage of various information about the observations. Firstly, these are text comments. Information about the place of work, the names of the companies that performed the measurements, etc. are entered here. In our case, this header is called by a request with the TextHeader key. Here's a shortened text heading:
Marmousi Model:
At the moment, the package also allows you to download data in parts from large files, making it possible to process files that can reach tens of gigabytes in size. The package also includes functions for exporting data to .segy and partial rewriting to the end of the file.
Separately, it is worth noting the functionality of the package when working with the complex structure of .segy files. Since it allows not only accessing keys and indices to individual traces, headers, but also changing them with subsequent writing to a file. Many of the technical details of implementing GeologyIO are beyond the scope of this article and probably deserve a separate description.
The ability to import seismic materials into Wolfram Mathematica allows you to use the built-in signal processing functionality for experimental data. Since each seismic track is a time series, one of the main tools for studying them is spectral analysis. Among the prerequisites for the analysis of the frequency composition of seismic data are, for example, the following:
The third paragraph is central in the context of this article. Below is a code snippet for calculating seismic traces in the case of a layer with a varying thickness - a wedge model. This model is traditionally studied in seismic exploration for the analysis of interference effects, when waves reflected from many layers overlap each other.
Outburst formation model:
The wave velocity inside the wedge is 4500 m / s, outside the wedge is 4000 m / s, and the density is assumed to be constant 2200 g / cm³. For such a model, we calculate the reflection coefficients and seismic traces.
Seismic tracks for the wedge model:
The sequence of seismic traces shown in this figure is called a seismic section. As you can see, its interpretation can be performed on an intuitive level, since the geometry of the reflected waves uniquely corresponds to the model that was set earlier. If you analyze the tracks in more detail, you can see that the tracks from the 1st to about the 30th do not differ - the reflection from the top of the formation and from the sole do not overlap. Starting from route 31, reflections begin to interfere. And, although, in the model, the reflection coefficients do not change horizontally - seismic tracks change their intensity with a change in the thickness of the reservoir.
Consider the amplitude of reflection from the upper boundary of the reservoir. Starting from the 60th route, the reflection intensity begins to increase and at the 70th route becomes maximum. This is how the interference of waves from the roof and sole for the strata manifests itself, leading in some cases to significant anomalies in the seismic record.
Graph of the amplitude of the reflected wave from the upper edge of the wedge
It is logical that when the signal is lower-frequency, then interference begins to appear at large thicknesses of the formation, and in the case of a high-frequency signal, interference occurs at smaller thicknesses. The following code fragment creates a signal with frequencies of 35 Hz, 55 Hz and 85 Hz.
A set of source signals with frequencies of 35 Hz, 55 Hz, 85 Hz
After calculating the seismic traces and plotting the amplitudes of the reflected wave, we can see that for different frequencies an anomaly is observed at different thicknesses of the formation.
Graphs of the amplitudes of the reflected wave from the upper edge of the wedge for different frequencies
The ability to draw conclusions about the thickness of the reservoir based on the results of seismic observations is extremely useful, because one of the main tasks in prospecting for an oil field is to evaluate the most promising points for laying a well (i.e., those areas where the reservoir has a large thickness). In addition, in the geological section, there can be found such objects that, due to their genesis, cause a sharp change in the thickness of the formation. This makes spectral analysis an effective tool for studying them. In the next part of the article, we will consider such gelogical objects in more detail.
The materials that are analyzed in the article were obtained on the territory of Western Siberia. The region, as everyone without exception probably knows, is the main oil-producing region of our country. Active development of metro-births began in the region in the 60s of the last century. The main method for finding oil fields is seismic exploration. It is interesting to consider satellite images of this territory. On a small scale, a huge number of swamps and lakes can be noted, by increasing the map, you can see cluster wells for drilling wells, and by increasing the map to the limit, you can also distinguish glades of profiles by which seismic observations were made.
Satellite image of Yandex maps - the area of the city of Noyabrsk:
A network of cluster sites at one of the fields:
The oil-bearing rocks of Western Siberia lie in a wide range of depths - from 1km to 5km. The main volume of rocks containing oil is formed in the Jurassic and Cretaceous times. The Jurassic period is probably known to many by the film of the same name. The climate of the Jurassic period was significantly different from the modern one. In Encyclopedia Britannica there is a series of paleocards that characterize every gelogic era.
Present time:
Jurassic period:
Please note that in the Jurassic, the territory of Western Siberia represented the sea coast (land crossed by rivers and shallow sea). Since the climate was comfortable, it can be assumed that a typical landscape of that time looked as follows:
Siberia of the Jurassic:
In this picture, not so much animals and birds are important for us, as the image of the river in the background. A river is the very geological object on which we stopped earlier. The fact is that the activity of rivers allows the accumulation of well-sorted sandstones, which then become a reservoir for oil. These reservoirs can have a bizarre, complex shape (like a river bed) and they have variable thickness - along the coast the thickness is small, and increases closer to the center of the channel or in sections of the meanders. So, the rivers formed in the Jurassic time are now at a depth of about three kilometers and are the object of the search for oil reservoirs.
We immediately make a reservation regarding the seismic materials shown in the article - due to the fact that the amount of data used for analysis is significant - only a fragment of the original set of seismic traces is placed in the text of the article. This will allow everyone to reproduce the above calculations.
When working with seismic data, geophysicists usually use specialized software (there are several industry leaders whose development is actively used, for example, Petrel or Paradigm), which allows you to analyze different types of data and has a convenient graphical interface. Despite all the convenience, these types of software also have their drawbacks - for example, the implementation of modern algorithms in stable versions takes a lot of time, and the possibilities of automation of calculations are usually limited. In such a situation, it becomes very convenient to use computer mathematics systems and high-level programming languages that allow you to use a wide algorithmic base and, at the same time, take on a lot of routine. Using this principle, work with seismic data in Wolfram Mathematica is built. It is inexpedient to write a rich functional of interactive work with data - it is more important to ensure loading from a generally accepted format, apply the desired algorithms to them and upload it back to an external format.
Following the proposed scheme, load the original seismic data and display it in Wolfram Mathematica :
The data downloaded and imported in this way are the tracks registered on a site measuring 10 by 5 kilometers. In the event that the data are obtained by the method of three-dimensional seismic exploration (wave registration is carried out not along separate geophysical profiles, but on the entire area at the same time), it becomes possible to obtain cubes of seismic data. These are three-dimensional objects, vertical and horizontal sections of which allow a detailed study of the geological environment. In the considered example, we are dealing with just three-dimensional data. We can get some information from the text heading, for example like this
One of the seismic sections:
Knowing how many traces are located in each direction of the studied area, you can form a three-dimensional data array and display it using the Image3D [] function
Three-dimensional image of a cube of seismic data. (Vertical axis - depth):
In the event that geological objects of interest create intense seismic anomalies, then visualization tools with transparency can be used. “Unimportant” sections of the recording can be made invisible, leaving only anomalies visible. In Wolfram Mathematica, this can be done using Opacity [] and Raster3D [] .
Image of a seismic data cube using the Opacity [] and Raster3D [] functions
As in the synthetic example, on the slices of the original cube, some geological boundaries (layers) with variable relief can be distinguished.
The main spectral analysis tool is the Fourier transform. With its help, you can evaluate the amplitude-frequency spectrum of each path or group of paths. However, after transferring data to the frequency domain, information is lost about at what times (read at what depths) the frequency changes. In order to be able to localize signal changes on the time (deep) axis, the window Fourier transform and wavelet decomposition are used. This article uses wavelet decomposition. Wavelet analysis technology began to be actively used in seismic exploration in the 90s. The advantage over the window Fourier transform is considered to be the best time resolution.
Using the following code fragment, you can decompose into individual components of one of the seismic traces:
Component decomposition
To assess how the reflection energy is distributed at different times of arrival of the waves, scalograms (analogue of the spectrogram) are used. As a rule, in practice, there is no need to analyze all components. Usually choose the low, medium and high frequency component.
Scalogram. Result of WaveletScalogram [] Function
Wolfram Language uses the ContinuousWaveletTransform [] function for wavelet transforms . And the application of this function to the entire set of traces is carried out with the use of the Table [] function . Here it is worth noting one of the strengths of Wolfram Mathematica is the ability to use ParallelTable [] parallelization . In the given example, parallelization is not necessary - the data volume is not large, but when working with experimental data sets containing hundreds of thousands of traces, this is a necessity.
After applying the ContinuousWaveletTransform [] function , new data arrays corresponding to the selected frequencies appear. In the above example, these are frequencies: 38Hz, 33Hz, 27Hz. The choice of frequencies is most often made on the basis of testing — they obtain effective maps for different frequency combinations and choose the most informative from the point of view of a geologist.
If you need to share the results with colleagues or provide them to the customer, you can use the SEGYExport [] function of GeologyIO
Having three such cubes available (low-frequency, mid-frequency and high-frequency components), as a rule, RGB mixing is used for joint visualization of data. Each component is assigned its own color - red, green, blue. In Wolfram Mathematica, this can be done using the ColorCombine [] function .
As a result, images are obtained that can be used to perform a geological interpretation. The meanders that hang on the section allow you to outline the paleorus, which is more likely to be reservoirs and contain oil reserves. Search and analysis of modern analogues of such a river system allows us to determine the most promising parts of the meanders. The beds themselves are characterized by thick layers of well-sorted sandstone and are a good reservoir for oil. Areas outside the "lace" anomalies are similar to modern floodplain deposits. Floodplain sediments are mainly represented by clayey rocks and drilling in these zones will be ineffective.
RGB slice of the data cube. In the center (a little to the left of the center) you can trace the meandering river.
RGB slice of the data cube. On the left side you can trace the meandering river.
In some cases, the quality of seismic data allows for significantly sharper images. It depends on the methods of field work, equipment, adopted by the noise reduction algorithm. In such cases, not only fragments of the river system are visible, but also entire extended paleorecs.
RGB mixing of the three components of the seismic data cube (horizontal slice). Depth is approximately 2 km.
Satellite image of the Volga river in the Saratov region
In Wolfram Mathematica, you can analyze seismic data and solve application problems associated with the search for minerals, and the GeologyIO package allows you to make this process more convenient. The structure of seismic data is such that the use of built-in methods for accelerating calculations ( ParallelTable [] , ParallelDo [], ...) is very efficient and allows you to process large amounts of data. To a large extent, this is facilitated by the data storage features of the GeologyIO package. By the way, the package can be used not only in the field of applied seismic exploration. Almost the same data types are used in georadar and seismology. If you have suggestions on how to improve the result, which Wolfram Mathematica signal analysis algorithms are applicable to such data, or if you have critical comments, leave comments.
Introduction
This note was written in the wake of the Wolfram Russian Technology Conference and contains a synopsis of the report I gave. The event took place in June in the city of St. Petersburg. Given the fact that I work in the quarter from the conference venue, I could not help but attend this event. In 2016 and 2017, I listened to the conference reports, and this year I made a presentation. Firstly, an interesting (it seems to me) topic has appeared that we are developing with Kirill Belov , and secondly, after a long study of the legislation of the Russian Federation regarding sanctions policy, two Wolfram Mathematica licenses have appeared in the enterprise where I work.
Before turning to the topic of my speech, I would like to note the good organization of the event. The conference page uses the image of the Kazan Cathedral. The Cathedral is one of the main attractions of St. Petersburg and it is very clearly visible from the hall in which the conference was held.
')
At the entrance to St. Petersburg State University of Economics, participants were met by student assistants - they did not allow them to get lost. During registration, small souvenirs were given out (a toy - flashing adhesions, a pen, stickers with Wolfram symbols). Lunch and coffee breaks were also included in the conference schedule. About delicious coffee and pies, I already noted on the wall of the group - well done chefs. With this introductory part, I would like to emphasize that the event itself, its format and venue already bring positive emotions.
The report, which was prepared by me and Kirill Belov, is called “Using Wolfram Mathematica to solve problems of applied geophysics. Spectral analysis of seismic data or "where the ancient rivers ran." The content of the report covers two parts: firstly, the use of algorithms available in Wolfram Mathematica for the analysis of geophysical data, and secondly, this is how to place geophysical data in Wolfram Mathematica.
Seismic exploration
First you need to do a little excursion into geophysics. Geophysics is a science that studies the physical properties of rocks. Well, since the rocks have different properties: electrical, magnetic, elastic, then there are appropriate methods of geophysics: electrical exploration, magnetic exploration, seismic exploration ... In the context of this article, we will only touch on seismic exploration. Seismic exploration is the main method for finding oil and gas. The method is based on the excitation of elastic vibrations and the subsequent recording of the response from rocks composing the study area. Excitation of vibrations is carried out on land (dynamite or non-explosive vibration sources of elastic vibrations) or at sea (air guns). Elastic vibrations propagate through the thickness of the rocks, being refracted and reflected at the boundaries of layers with different properties. The reflected waves return to the surface and are recorded by land-based geophones (usually electrodynamic devices based on the movement of a magnet suspended in a coil) or hydrophones in the sea (based on the piezoelectric effect). By the time of arrival of the waves, one can judge the depths of the geological layers.
Seismic vessel tows equipment:
Air gun excites elastic vibrations:
Waves pass through the rock mass and are recorded by hydrophones:
Geophysical exploration research vessel “Ivan Gubkin” on the pier near the Blagoveshchensky bridge in St. Petersburg:
Seismic signal model
Rocks have different physical properties. For seismic exploration, elastic properties — the propagation velocity of elastic vibrations and density — are of primary importance. If two layers have the same or similar properties, the wave will not “notice” the boundary between them. If the wave velocities in the layers differ, then reflection will occur at the boundary of the layers. The greater the difference in properties, the more intense the reflection. Its intensity will be determined by the reflection coefficient (rc):
where ρ is the rock density, ν is the wave velocity, 1 and 2 denote the upper and lower layers.
One of the simplest and most commonly used seismic signal models is the convolution model, when a registered seismic trace is presented as a result of convolution of a sequence of reflection coefficients with a probe pulse:
where s (t) is the seismic track, i.e. all that a hydrophone or a geophone recorded over a fixed recording time, w (t) is the signal that the airgun generates, n (t) is random noise.
Let us calculate, for example, a synthetic seismic survey. We will use the Ricker pulse widely used in seismic exploration as the initial signal.
length=0.050; (*Signal lenght*) dt=0.001;(*Sample rate of signal*) t=Range[-length/2,(length)/2,dt];(*Signal time*) f=35;(*Central frequency*) wavelet=(1.0-2.0*(Pi^2)*(f^2)*(t^2))*Exp[-(Pi^2)*(f^2)*(t^2)]; ListLinePlot[wavelet, Frame->True,PlotRange->Full,Filling->Axis,PlotStyle->Black, PlotLabel->Style["Initial wavelet",Black,20], LabelStyle->Directive[Black,Italic], FillingStyle->{White,Black},ImageSize->Large,InterpolationOrder->2] Source seismic pulse
We set two boundaries at depths of 300 ms and 600 ms, and the reflection coefficients will be random numbers.
rcExample=ConstantArray[0,1000]; rcExample[[300]]=RandomReal[{-1,0}]; rcExample[[600]]=RandomReal[{0,1}]; ListPlot[rcExample,Filling->0,Frame->True,Axes->False,PlotStyle->Black, PlotLabel->Style["Reflection Coefficients",Black,20], LabelStyle->Directive[Black,Italic]] The sequence of reflection coefficients:
Calculate and display the seismic track. Since the reflection coefficients have different signs, we obtain two alternating reflections on the seismic path.
traceExamle=ListConvolve[wavelet[[1;;;;1]],rcExample]; ListPlot[traceExamle, PlotStyle->Black,Filling->0,Frame->True,Axes->False, PlotLabel->Style["Seismic trace",Black,20], LabelStyle->Directive[Black,Italic]] Simulated track:
For this example, a reservation must be made - in reality, the depth of the layers is determined of course in meters, and the calculation of the seismic path occurs for the time domain. It would be more correct to set the depths in meters and calculate the arrival times knowing the velocity in the strata. In this case, I immediately set the layers on the time axis.
If we talk about field studies, as a result of such observations, a huge number of such time series (seismic traces) are recorded. For example, when exploring a site 25 km long and 15 km wide, where as a result of work each track characterizes a cell measuring 25x25 meters (this cell is called a bin), the final data array will contain 600,000 tracks. With a time sampling step of 1 ms, a recording time of 5 seconds, the final data file will be more than 11 GB, and the amount of the original “raw” material can be hundreds of gigabytes.
How to work with them in Wolfram Mathematica ?
GeologyIO Package
The development of the package began with a question on the VK wall of the Russian-speaking support group. Thanks to community responses, a solution was found very quickly. And as a result, it grew into a serious development. The corresponding post on the wall of Wolfram Comunity was even marked by moderators. Currently, the package supports working with the following types of data that are actively used in the geological industry:
- import map data of ZMAP and IRAP format
- import of measurements in LAS format wells
- input and output of seismic files in SEGY format
To install the package, you must follow the instructions on the download page of the assembled package, i.e. execute the following code on any Mathematica notebook :
If[PacletInformation["GeologyIO"] === {}, PacletInstall[URLDownload[ "https://wolfr.am/FiQ5oFih", FileNameJoin[{CreateDirectory[], "GeologyIO-0.2.2.paclet"}] ]]] After which the package will be installed in the default folder, the path to which can be obtained as follows:
FileNameJoin[{$UserBasePacletsDirectory, "Repository"}] For example, we demonstrate the main features of the package. The call is traditionally made for packages in Wolfram Language:
Get["GeologyIO`"] The package is being developed using Wolfram Workbench . This allows you to accompany the main functionality of the package with documentation, which in presentation format does not differ from the documentation of Wolfram Mathematica itself and provide the package with test files for the first acquaintance.
Such a file, in particular, is the Marmousi.segy file, a synthetic model of a geological section developed by the French Petroleum Institute. Using this model, developers test their own algorithms for modeling the wave field, data processing, inversion of seismic traces, etc. The Marmousi model itself is stored in the repository, from where the package itself was downloaded. In order to get the file, execute the following code:
If[Not[FileExistsQ["Marmousi.segy"]], URLDownload["https://wolfr.am/FiQGh7rk", "Marmousi.segy"];] marmousi = SEGYImport["Marmousi.segy"] The result of the import is a SEGYData object:
The SEGY format involves the storage of various information about the observations. Firstly, these are text comments. Information about the place of work, the names of the companies that performed the measurements, etc. are entered here. In our case, this header is called by a request with the TextHeader key. Here's a shortened text heading:
Short[marmousi["TextHeader"]] “The Marmousi data set was generated at the Institute ... nimum velocity of 1500 m / s and a maximum of 5500 m / s)”The actual geological model can be displayed by contacting the seismic traces using the “traces” key (one of the features of the package is the case-insensitivity of the keys):
ArrayPlot[Transpose[marmousi["traces"]], PlotTheme -> "Detailed"] Marmousi Model:
At the moment, the package also allows you to download data in parts from large files, making it possible to process files that can reach tens of gigabytes in size. The package also includes functions for exporting data to .segy and partial rewriting to the end of the file.
Separately, it is worth noting the functionality of the package when working with the complex structure of .segy files. Since it allows not only accessing keys and indices to individual traces, headers, but also changing them with subsequent writing to a file. Many of the technical details of implementing GeologyIO are beyond the scope of this article and probably deserve a separate description.
The relevance of spectral analysis in seismic exploration
The ability to import seismic materials into Wolfram Mathematica allows you to use the built-in signal processing functionality for experimental data. Since each seismic track is a time series, one of the main tools for studying them is spectral analysis. Among the prerequisites for the analysis of the frequency composition of seismic data are, for example, the following:
- Different types of waves are characterized by different frequency composition. This allows us to select useful waves and suppress interference waves.
- Rock properties such as porosity and saturation can affect the frequency composition. This allows you to select rocks with better properties.
- Layers with different thicknesses cause anomalies in different frequency ranges.
The third paragraph is central in the context of this article. Below is a code snippet for calculating seismic traces in the case of a layer with a varying thickness - a wedge model. This model is traditionally studied in seismic exploration for the analysis of interference effects, when waves reflected from many layers overlap each other.
nx=200;(* Number of grid points in X direction*) ny=200;(* Number of grid points in Y direction*) T=2;(*Total propagation time*) (*Velocity and density*) modellv=Table[4000,{i,1,ny},{j,1,nx}];(* P-wave velocity in m/s*) rho=Table[2200,{i,1,ny},{j,1,nx}];(* Density in g/cm^3, used constant density*) Table[modellv[[150-Round[i*0.5];;,i]]=4500;,{i,1,200}]; Table[modellv[[;;70,i]]=4500;,{i,1,200}]; (*Plotting model*) MatrixPlot[modellv,PlotLabel->Style["Model of layer",Black,20], LabelStyle->Directive[Black,Italic]] Outburst formation model:
The wave velocity inside the wedge is 4500 m / s, outside the wedge is 4000 m / s, and the density is assumed to be constant 2200 g / cm³. For such a model, we calculate the reflection coefficients and seismic traces.
rc=Table[N[(modellv[[All,i]]-PadLeft[modellv[[All,i]],201,4000][[1;;200]])/(modellv[[All,i]]+PadLeft[modellv[[All,i]],201,4500][[1;;200]])],{i,1,200}]; traces=Table[ListConvolve[wavelet[[1;;;;1]],rc[[i]]],{i,1,200}]; starttrace=10; endtrace=200; steptrace=10; trasenum=Range[starttrace,endtrace,steptrace]; traserenum=Range[Length@trasenum]; tracedist=0.5; Rotate[Show[ Reverse[Table[ ListLinePlot[traces[[trasenum[[i]]]]*50+trasenum[[i]]*tracedist,Filling->{1->{trasenum[[i]]*tracedist,{RGBColor[0.97,0.93,0.68],Black}}},PlotStyle->Directive[Gray,Thin],PlotRange->Full,InterpolationOrder->2,Axes->False,Background->RGBColor[0.97,0.93,0.68]], {i,1,Length@trasenum}]],ListLinePlot[Transpose[{ConstantArray[45,80],Range[80]}],PlotStyle->Red],PlotRange->All,Frame->True],270Degree] Seismic tracks for the wedge model:
The sequence of seismic traces shown in this figure is called a seismic section. As you can see, its interpretation can be performed on an intuitive level, since the geometry of the reflected waves uniquely corresponds to the model that was set earlier. If you analyze the tracks in more detail, you can see that the tracks from the 1st to about the 30th do not differ - the reflection from the top of the formation and from the sole do not overlap. Starting from route 31, reflections begin to interfere. And, although, in the model, the reflection coefficients do not change horizontally - seismic tracks change their intensity with a change in the thickness of the reservoir.
Consider the amplitude of reflection from the upper boundary of the reservoir. Starting from the 60th route, the reflection intensity begins to increase and at the 70th route becomes maximum. This is how the interference of waves from the roof and sole for the strata manifests itself, leading in some cases to significant anomalies in the seismic record.
ListLinePlot[GaussianFilter[Abs[traces[[All,46]]],3][[;;;;2]], InterpolationOrder->2,Frame->True,PlotStyle->Black, PlotLabel->Style["Amplitude of reflection",Black,20], LabelStyle->Directive[Black,Italic], PlotRange->All] Graph of the amplitude of the reflected wave from the upper edge of the wedge
It is logical that when the signal is lower-frequency, then interference begins to appear at large thicknesses of the formation, and in the case of a high-frequency signal, interference occurs at smaller thicknesses. The following code fragment creates a signal with frequencies of 35 Hz, 55 Hz and 85 Hz.
waveletSet=Table[(1.0-2.0*(Pi^2)*(f^2)*(t^2))*Exp[-(Pi^2)*(f^2)*(t^2)], {f,{35,55,85}}]; ListLinePlot[waveletSet,PlotRange->Full,PlotStyle->Black,Frame->True, PlotLabel->Style["Set of wavelets",Black,20], LabelStyle->Directive[Black,Italic], ImageSize->Large,InterpolationOrder->2] A set of source signals with frequencies of 35 Hz, 55 Hz, 85 Hz
After calculating the seismic traces and plotting the amplitudes of the reflected wave, we can see that for different frequencies an anomaly is observed at different thicknesses of the formation.
tracesSet=Table[ListConvolve[waveletSet[[j]][[1;;;;1]],rc[[i]]],{j,1,3},{i,1,200}]; lowFreq=ListLinePlot[GaussianFilter[Abs[tracesSet[[1]][[All,46]]],3][[;;;;2]],InterpolationOrder->2,PlotStyle->Black,PlotRange->All]; medFreq=ListLinePlot[GaussianFilter[Abs[tracesSet[[2]][[All,46]]],3][[;;;;2]],InterpolationOrder->2,PlotStyle->Black,PlotRange->All]; highFreq=ListLinePlot[GaussianFilter[Abs[tracesSet[[3]][[All,46]]],3][[;;;;2]],InterpolationOrder->2,PlotStyle->Black,PlotRange->All]; Show[lowFreq,medFreq,highFreq,PlotRange->{{0,100},All}, PlotLabel->Style["Amplitudes of reflection",Black,20], LabelStyle->Directive[Black,Italic], Frame->True] Graphs of the amplitudes of the reflected wave from the upper edge of the wedge for different frequencies
The ability to draw conclusions about the thickness of the reservoir based on the results of seismic observations is extremely useful, because one of the main tasks in prospecting for an oil field is to evaluate the most promising points for laying a well (i.e., those areas where the reservoir has a large thickness). In addition, in the geological section, there can be found such objects that, due to their genesis, cause a sharp change in the thickness of the formation. This makes spectral analysis an effective tool for studying them. In the next part of the article, we will consider such gelogical objects in more detail.
Experimental data. Where are they received and what to look for in them?
The materials that are analyzed in the article were obtained on the territory of Western Siberia. The region, as everyone without exception probably knows, is the main oil-producing region of our country. Active development of metro-births began in the region in the 60s of the last century. The main method for finding oil fields is seismic exploration. It is interesting to consider satellite images of this territory. On a small scale, a huge number of swamps and lakes can be noted, by increasing the map, you can see cluster wells for drilling wells, and by increasing the map to the limit, you can also distinguish glades of profiles by which seismic observations were made.
Satellite image of Yandex maps - the area of the city of Noyabrsk:
A network of cluster sites at one of the fields:
The oil-bearing rocks of Western Siberia lie in a wide range of depths - from 1km to 5km. The main volume of rocks containing oil is formed in the Jurassic and Cretaceous times. The Jurassic period is probably known to many by the film of the same name. The climate of the Jurassic period was significantly different from the modern one. In Encyclopedia Britannica there is a series of paleocards that characterize every gelogic era.
Present time:
Jurassic period:
Please note that in the Jurassic, the territory of Western Siberia represented the sea coast (land crossed by rivers and shallow sea). Since the climate was comfortable, it can be assumed that a typical landscape of that time looked as follows:
Siberia of the Jurassic:
In this picture, not so much animals and birds are important for us, as the image of the river in the background. A river is the very geological object on which we stopped earlier. The fact is that the activity of rivers allows the accumulation of well-sorted sandstones, which then become a reservoir for oil. These reservoirs can have a bizarre, complex shape (like a river bed) and they have variable thickness - along the coast the thickness is small, and increases closer to the center of the channel or in sections of the meanders. So, the rivers formed in the Jurassic time are now at a depth of about three kilometers and are the object of the search for oil reservoirs.
Experimental data. Processing and visualization
We immediately make a reservation regarding the seismic materials shown in the article - due to the fact that the amount of data used for analysis is significant - only a fragment of the original set of seismic traces is placed in the text of the article. This will allow everyone to reproduce the above calculations.
When working with seismic data, geophysicists usually use specialized software (there are several industry leaders whose development is actively used, for example, Petrel or Paradigm), which allows you to analyze different types of data and has a convenient graphical interface. Despite all the convenience, these types of software also have their drawbacks - for example, the implementation of modern algorithms in stable versions takes a lot of time, and the possibilities of automation of calculations are usually limited. In such a situation, it becomes very convenient to use computer mathematics systems and high-level programming languages that allow you to use a wide algorithmic base and, at the same time, take on a lot of routine. Using this principle, work with seismic data in Wolfram Mathematica is built. It is inexpedient to write a rich functional of interactive work with data - it is more important to ensure loading from a generally accepted format, apply the desired algorithms to them and upload it back to an external format.
Following the proposed scheme, load the original seismic data and display it in Wolfram Mathematica :
Get["GeologyIO`"] seismic3DZipPath = "seismic3D.zip"; seismic3DSEGYPath = "seismic3D.sgy"; If[FileExistsQ[seismic3DZipPath], DeleteFile[seismic3DZipPath]]; If[FileExistsQ[seismic3DSEGYPath], DeleteFile[seismic3DSEGYPath]]; URLDownload["https://wolfr.am/FiQIuZuH", seismic3DZipPath]; ExtractArchive[seismic3DZipPath]; seismic3DSEGY = SEGYImport[seismic3DSEGYPath] The data downloaded and imported in this way are the tracks registered on a site measuring 10 by 5 kilometers. In the event that the data are obtained by the method of three-dimensional seismic exploration (wave registration is carried out not along separate geophysical profiles, but on the entire area at the same time), it becomes possible to obtain cubes of seismic data. These are three-dimensional objects, vertical and horizontal sections of which allow a detailed study of the geological environment. In the considered example, we are dealing with just three-dimensional data. We can get some information from the text heading, for example like this
StringPartition[seismic3DSEGY["textheader"], 80] // TableForm C 1 THIS IS DEMO FILE FOR GEOLOGYIO PACKAGE TESTThis data set will be enough for us to demonstrate the main stages of data analysis. The tracks in the file are recorded sequentially and each of them looks something like the following figure — this is the distribution of the amplitudes of the reflected waves along the vertical axis (depth axis).
C 2
C 3
C 4
C 5 DATE USER NAME: WOLFRAM USER
C 6 SURVEY NAME: SOMEWHERE IN SIBERIA
C 7 FILE TYPE 3D SEISMIC VOLUME
C 8
C 9
C10 Z RANGE: FIRST 2200M LAST 2400M
ListLinePlot[seismic3DSEGY["traces"][[100]], InterpolationOrder -> 2, PlotStyle -> Black, PlotLabel -> Style["Seismic trace", Black, 20], LabelStyle -> Directive[Black, Italic], PlotRange -> All, Frame -> True, ImageSize -> 1200, AspectRatio -> 1/5] One of the seismic sections:
Knowing how many traces are located in each direction of the studied area, you can form a three-dimensional data array and display it using the Image3D [] function
traces=seismic3DSEGY["traces"]; startIL=1050;EndIL=2000;stepIL=2; (* *) startXL=1165;EndXL=1615;stepXL=2; (* Y *) numIL=(EndIL-startIL)/stepIL+1; (* *) numXL=(EndXL-startXL)/stepIL+1; (* Y*) Image3D[ArrayReshape[Abs[traces/Max[Abs[traces[[All,1;;;;4]]]]],{numIL,numXL,101}],ViewPoint->{-1, 0, 0},Background->RGBColor[0,0,0]] Three-dimensional image of a cube of seismic data. (Vertical axis - depth):
In the event that geological objects of interest create intense seismic anomalies, then visualization tools with transparency can be used. “Unimportant” sections of the recording can be made invisible, leaving only anomalies visible. In Wolfram Mathematica, this can be done using Opacity [] and Raster3D [] .
data = ArrayReshape[Abs[traces/Max[Abs[traces[[All,1;;;;4]]]]],{numIL,numXL,101}]; Graphics3D[{Opacity[0.1], Raster3D[data, ColorFunction->"RainbowOpacity"]}, Boxed->False, SphericalRegion->True, ImageSize->840, Background->None] Image of a seismic data cube using the Opacity [] and Raster3D [] functions
As in the synthetic example, on the slices of the original cube, some geological boundaries (layers) with variable relief can be distinguished.
The main spectral analysis tool is the Fourier transform. With its help, you can evaluate the amplitude-frequency spectrum of each path or group of paths. However, after transferring data to the frequency domain, information is lost about at what times (read at what depths) the frequency changes. In order to be able to localize signal changes on the time (deep) axis, the window Fourier transform and wavelet decomposition are used. This article uses wavelet decomposition. Wavelet analysis technology began to be actively used in seismic exploration in the 90s. The advantage over the window Fourier transform is considered to be the best time resolution.
Using the following code fragment, you can decompose into individual components of one of the seismic traces:
cwd=ContinuousWaveletTransform[seismicSection["traces"][[100]]] Show[ ListLinePlot[Re[cwd[[1]]],PlotRange->All], ListLinePlot[seismicSection["traces"][[100]], PlotStyle->Black,PlotRange->All],ImageSize->{1500,500},AspectRatio->Full, PlotLabel->Style["Wavelet decomposition",Black,32], LabelStyle->Directive[Black,Italic], PlotRange->All, Frame->True] Component decomposition
To assess how the reflection energy is distributed at different times of arrival of the waves, scalograms (analogue of the spectrogram) are used. As a rule, in practice, there is no need to analyze all components. Usually choose the low, medium and high frequency component.
freq=(500/(#*contWD["Wavelet"]["FourierFactor"]))&/@(Thread[{Range[contWD["Octaves"]],1}]/.contWD["Scales"])//Round; ticks=Transpose[{Range[Length[freq]],freq}]; WaveletScalogram[contWD,Frame->True,FrameTicks->{{ticks,Automatic},Automatic},FrameTicksStyle->Directive[Orange,12], FrameLabel->{"Time","Frequency(Hz)"},LabelStyle->Directive[Black,Bold,14], ColorFunction->"RustTones",ImageSize->Large] Scalogram. Result of WaveletScalogram [] Function
Wolfram Language uses the ContinuousWaveletTransform [] function for wavelet transforms . And the application of this function to the entire set of traces is carried out with the use of the Table [] function . Here it is worth noting one of the strengths of Wolfram Mathematica is the ability to use ParallelTable [] parallelization . In the given example, parallelization is not necessary - the data volume is not large, but when working with experimental data sets containing hundreds of thousands of traces, this is a necessity.
tracesCWD=Table[Map[Hilbert[#,0]&,Re[ContinuousWaveletTransform[traces[[i]]][[1]]][[{13,15,18}]]],{i,1,Length@traces}]; After applying the ContinuousWaveletTransform [] function , new data arrays corresponding to the selected frequencies appear. In the above example, these are frequencies: 38Hz, 33Hz, 27Hz. The choice of frequencies is most often made on the basis of testing — they obtain effective maps for different frequency combinations and choose the most informative from the point of view of a geologist.
If you need to share the results with colleagues or provide them to the customer, you can use the SEGYExport [] function of GeologyIO
outputdata=seismic3DSEGY; outputdata["traces",1;;-1]=tracesCWD[[All,3]]; outputdata["textheader"]="Wavelet Decomposition Result"; outputdata["binaryheader","NumberDataTraces"]=Length[tracesCWD[[All,3]]]; SEGYExport["D:\\result.segy",outputdata]; Having three such cubes available (low-frequency, mid-frequency and high-frequency components), as a rule, RGB mixing is used for joint visualization of data. Each component is assigned its own color - red, green, blue. In Wolfram Mathematica, this can be done using the ColorCombine [] function .
As a result, images are obtained that can be used to perform a geological interpretation. The meanders that hang on the section allow you to outline the paleorus, which is more likely to be reservoirs and contain oil reserves. Search and analysis of modern analogues of such a river system allows us to determine the most promising parts of the meanders. The beds themselves are characterized by thick layers of well-sorted sandstone and are a good reservoir for oil. Areas outside the "lace" anomalies are similar to modern floodplain deposits. Floodplain sediments are mainly represented by clayey rocks and drilling in these zones will be ineffective.
RGB slice of the data cube. In the center (a little to the left of the center) you can trace the meandering river.
RGB slice of the data cube. On the left side you can trace the meandering river.
In some cases, the quality of seismic data allows for significantly sharper images. It depends on the methods of field work, equipment, adopted by the noise reduction algorithm. In such cases, not only fragments of the river system are visible, but also entire extended paleorecs.
RGB mixing of the three components of the seismic data cube (horizontal slice). Depth is approximately 2 km.
Satellite image of the Volga river in the Saratov region
Conclusion
In Wolfram Mathematica, you can analyze seismic data and solve application problems associated with the search for minerals, and the GeologyIO package allows you to make this process more convenient. The structure of seismic data is such that the use of built-in methods for accelerating calculations ( ParallelTable [] , ParallelDo [], ...) is very efficient and allows you to process large amounts of data. To a large extent, this is facilitated by the data storage features of the GeologyIO package. By the way, the package can be used not only in the field of applied seismic exploration. Almost the same data types are used in georadar and seismology. If you have suggestions on how to improve the result, which Wolfram Mathematica signal analysis algorithms are applicable to such data, or if you have critical comments, leave comments.
Source: https://habr.com/ru/post/461241/
All Articles