Gracias al autor del blog. por su informe
introducción
Esta nota fue escrita a raíz de la conferencia. y contiene un resumen del informe que di. El evento tuvo lugar en junio en San Petersburgo. Considerando que trabajo a una cuadra del lugar de la conferencia, no pude evitar asistir a este evento. En 2016 y 2017 escuché informes de conferencias y este año hice una presentación. En primer lugar, ha aparecido un tema interesante (me parece) que estamos desarrollando con , y en segundo lugar, después de un largo estudio de la legislación de la Federación de Rusia sobre la política de sanciones, en la empresa donde trabajo aparecieron hasta dos licencias. .
Antes de pasar al tema de mi intervención, quisiera destacar la buena organización del evento. La página de visitas de la conferencia utiliza una imagen de la Catedral de Kazán. La catedral es una de las principales atracciones de San Petersburgo y es muy visible desde la sala en la que tuvo lugar la conferencia.
A la entrada de la Universidad Estatal de Economía de San Petersburgo, los participantes fueron recibidos por asistentes de entre los estudiantes que no les permitieron perderse. Durante el registro se repartieron pequeños souvenirs (un juguete: una púa parpadeante, un bolígrafo, pegatinas con los símbolos de Wolfram). El almuerzo y las pausas para el café también se incluyeron en el programa de la conferencia. Ya mencioné los deliciosos cafés y pasteles en la pared del grupo: los chefs son geniales. Con esta parte introductoria, me gustaría enfatizar que el evento en sí, su formato y ubicación ya están generando emociones positivas.
El informe que preparamos Kirill Belov y yo se llama “Uso de Wolfram Mathematica para resolver problemas en geofísica aplicada”. Análisis espectral de datos sísmicos o “dónde corrían ríos antiguos”. El contenido del informe cubre dos partes: en primer lugar, el uso de algoritmos disponibles en para analizar datos geofísicos y, en segundo lugar, cómo colocar datos geofísicos en Wolfram Mathematica.
Exploración sísmica
Primero necesitas hacer una breve excursión a la geofísica. La geofísica es la ciencia que estudia las propiedades físicas de las rocas. Bueno, dado que las rocas tienen diferentes propiedades: eléctrica, magnética, elástica, existen métodos geofísicos correspondientes: prospección eléctrica, prospección magnética, prospección sísmica... En el contexto de este artículo, solo discutiremos la prospección sísmica con más detalle. La exploración sísmica es el principal método de búsqueda de petróleo y gas. El método se basa en la excitación de vibraciones elásticas y el posterior registro de la respuesta de las rocas que componen el área de estudio. Las vibraciones se excitan en tierra (con dinamita o fuentes de vibración elásticas no explosivas) o en el mar (con pistolas de aire comprimido). Las vibraciones elásticas se propagan a través del macizo rocoso, siendo refractadas y reflejadas en los límites de capas con diferentes propiedades. Las ondas reflejadas regresan a la superficie y son registradas por geófonos en tierra (normalmente dispositivos electrodinámicos basados en el movimiento de un imán suspendido en una bobina) o hidrófonos en el mar (basados en el efecto piezoeléctrico). Por el momento de la llegada de las olas, se puede juzgar la profundidad de las capas geológicas.
Equipo de remolque de buques sísmicos.
La pistola de aire excita vibraciones elásticas.
Las ondas atraviesan el macizo rocoso y son registradas por hidrófonos.
El buque de investigación geofísica "Ivan Gubkin" en el muelle cerca del puente Blagoveshchensky en San Petersburgo
Modelo de señal sísmica
Las rocas tienen diferentes propiedades físicas. Para la exploración sísmica, las propiedades elásticas son principalmente importantes: la velocidad de propagación de las vibraciones elásticas y la densidad. Si dos capas tienen propiedades iguales o similares, entonces la onda "no notará" el límite entre ellas. Si las velocidades de las ondas en las capas difieren, entonces se producirá una reflexión en el límite de las capas. Cuanto mayor es la diferencia de propiedades, más intensa es la reflexión. Su intensidad vendrá determinada por el coeficiente de reflectancia (rc):
donde ρ es la densidad de la roca, ν es la velocidad de la onda, 1 y 2 indican las capas superior e inferior.
Uno de los modelos de señales sísmicas más simples y más utilizados es el modelo de convolución, cuando la traza sísmica registrada se representa como resultado de la convolución de una secuencia de coeficientes de reflexión con un pulso de sondeo:
donde s(t) — rastro sísmico, es decir todo lo que fue registrado por un hidrófono o geófono durante un tiempo de grabación fijo, peso(t) - la señal generada por la pistola de aire comprimido, Nuevo Testamento) - ruido aleatorio.
Calculemos como ejemplo una traza sísmica sintética. Utilizaremos como señal inicial el pulso de Ricker, muy utilizado en exploración sísmica.
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]Impulso sísmico inicial
Estableceremos dos límites a profundidades de 300 ms y 600 ms, y los coeficientes de reflexión serán números aleatorios.
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]]Secuencia de coeficientes de reflexión.
Calculemos y visualicemos la traza sísmica. Dado que los coeficientes de reflexión tienen signos diferentes, obtenemos dos reflexiones alternas sobre la traza sísmica.
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]]Pista simulada
Para este ejemplo, es necesario hacer una reserva: en realidad, la profundidad de las capas se determina, por supuesto, en metros, y el cálculo de la traza sísmica se realiza para el dominio del tiempo. Sería más correcto fijar las profundidades en metros y calcular los tiempos de llegada conociendo las velocidades en las capas. En este caso, inmediatamente coloqué las capas en el eje del tiempo.
Si hablamos de investigación de campo, como resultado de tales observaciones se registra una gran cantidad de series temporales similares (trazas sísmicas). Por ejemplo, al estudiar un sitio de 25 km de largo y 15 km de ancho, donde, como resultado del trabajo, cada rastro caracteriza una celda de 25x25 metros (dicha celda se llama contenedor), la matriz de datos final contendrá 600000 rastros. Con un tiempo de muestreo de 1 ms y un tiempo de grabación de 5 segundos, el archivo de datos final tendrá más de 11 GB y el volumen del material "en bruto" original puede ser de cientos de gigabytes.
Cómo trabajar con ellos ?
Paquete
El desarrollo del paquete comenzó. en el muro VK del grupo de apoyo de habla rusa. Gracias a las respuestas de la comunidad, se encontró una solución muy rápidamente. Y como resultado, se convirtió en un desarrollo serio. Correspondiente Incluso fue marcado por los moderadores. Actualmente, el paquete admite trabajar con los siguientes tipos de datos que se utilizan activamente en la industria geológica:
- importación de datos de mapas en formatos ZMAP e IRAP
- importación de mediciones en pozos formato LAS
- Formato de entrada y salida de archivos sísmicos.
Para instalar el paquete, debe seguir las instrucciones en la página de descarga del paquete ensamblado, es decir ejecuta el siguiente código en cualquier :
If[PacletInformation["GeologyIO"] === {}, PacletInstall[URLDownload[
"https://wolfr.am/FiQ5oFih",
FileNameJoin[{CreateDirectory[], "GeologyIO-0.2.2.paclet"}]
]]]Después de lo cual el paquete se instalará en la carpeta predeterminada, cuya ruta se puede obtener de la siguiente manera:
FileNameJoin[{$UserBasePacletsDirectory, "Repository"}]Como ejemplo, demostraremos las principales capacidades del paquete. La llamada se realiza tradicionalmente para paquetes en Wolfram Language:
Get["GeologyIO`"]El paquete se desarrolla utilizando . Esto le permite acompañar la funcionalidad principal del paquete con documentación, que en términos de formato de presentación no difiere de la documentación de Wolfram Mathematica, y proporcionar al paquete archivos de prueba para su primer conocimiento.
Un archivo de este tipo, en particular, es el archivo "Marmousi.segy": se trata de un modelo sintético de una sección geológica desarrollado por el Instituto Francés del Petróleo. Utilizando este modelo, los desarrolladores prueban sus propios algoritmos para modelado de campos de ondas, procesamiento de datos, inversión de trazas sísmicas, etc. El modelo Marmousi en sí se almacena en el repositorio desde donde se descargó el paquete. Para obtener el archivo, ejecute el siguiente código:
If[Not[FileExistsQ["Marmousi.segy"]],
URLDownload["https://wolfr.am/FiQGh7rk", "Marmousi.segy"];]
marmousi = SEGYImport["Marmousi.segy"]Resultado de la importación: objeto SEGYData
El formato SEGY implica almacenar información diversa sobre las observaciones. En primer lugar, estos son comentarios de texto. Esto incluye información sobre la ubicación de la obra, los nombres de las empresas que realizaron las mediciones, etc. En nuestro caso, este encabezado se llama mediante una solicitud con la clave TextHeader. Aquí hay un encabezado de texto abreviado:
Short[marmousi["TextHeader"]]“El conjunto de datos de Marmousi fue generado en el Instituto...velocidad mínima de 1500 m/s y máxima de 5500 m/s)”
Puede mostrar el modelo geológico real accediendo a los rastros sísmicos usando la tecla "traces" (una de las características del paquete es que las claves no distinguen entre mayúsculas y minúsculas):
ArrayPlot[Transpose[marmousi["traces"]], PlotTheme -> "Detailed"]Modelo Marmusi
Actualmente, el paquete también permite cargar datos en partes desde archivos grandes, lo que permite procesar archivos cuyo tamaño puede alcanzar decenas de gigabytes. Las funciones del paquete también incluyen funciones para exportar datos a .segy y agregarlos parcialmente al final del archivo.
Por otra parte, cabe destacar la funcionalidad del paquete cuando se trabaja con la compleja estructura de archivos .segy. Ya que le permite no solo acceder a rastros y encabezados individuales usando claves e índices, sino también cambiarlos y luego escribirlos en un archivo. Muchos de los detalles técnicos de la implementación de GeologyIO están fuera del alcance de este artículo y probablemente merezcan una descripción por separado.
Relevancia del análisis espectral en la exploración sísmica.
La capacidad de importar datos sísmicos a Wolfram Mathematica le permite utilizar la funcionalidad de procesamiento de señales incorporada para datos experimentales. Dado que cada traza sísmica representa una serie temporal, una de las principales herramientas para estudiarlas es el análisis espectral. Entre los requisitos previos para analizar la composición frecuencial de los datos sísmicos podemos nombrar, por ejemplo, los siguientes:
- Los diferentes tipos de ondas se caracterizan por una composición de frecuencia diferente. Esto le permite resaltar ondas útiles y suprimir ondas de interferencia.
- Las propiedades de la roca, como la porosidad y la saturación, pueden afectar la composición de frecuencia. Esto permite identificar rocas con las mejores propiedades.
- Las capas con diferentes espesores provocan anomalías en diferentes rangos de frecuencia.
El tercer punto es el principal en el contexto de este artículo. A continuación se muestra un fragmento de código para calcular las trazas sísmicas en el caso de una capa con espesor variable: un modelo de cuña. Este modelo se estudia tradicionalmente en la exploración sísmica para analizar los efectos de interferencia cuando las ondas reflejadas de muchas capas se superponen entre sí.
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]]Modelo de formación pinch-out
La velocidad de la onda dentro de la cuña es de 4500 m/s, fuera de la cuña de 4000 m/s y se supone que la densidad es constante de 2200 g/cm³. Para tal modelo, calculamos coeficientes de reflexión y trazas sísmicas.
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]Trazas sísmicas para el modelo de cuña.
La secuencia de trazas sísmicas que se muestra en esta figura se denomina sección sísmica. Como puede ver, su interpretación también se puede realizar a nivel intuitivo, ya que la geometría de las ondas reflejadas corresponde claramente al modelo especificado anteriormente. Si analiza las huellas con más detalle, notará que las huellas del 1 al 30 aproximadamente no difieren: el reflejo desde el techo de la formación y desde abajo no se superponen entre sí. A partir del trazo 31, los reflejos comienzan a interferir. Y, aunque en el modelo los coeficientes de reflexión no cambian horizontalmente, las huellas sísmicas cambian de intensidad a medida que cambia el espesor de la formación.
Consideremos la amplitud de reflexión desde el límite superior de la formación. A partir de la ruta 60, la intensidad del reflejo comienza a aumentar y en la ruta 70 se vuelve máxima. Así se manifiesta la interferencia de ondas provenientes del techo y del fondo de las capas, provocando en algunos casos anomalías importantes en el registro sísmico.
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]Gráfico de la amplitud de la onda reflejada desde el borde superior de la cuña.
Es lógico que cuando la señal es de baja frecuencia, la interferencia comience a aparecer en espesores de formación grandes, y en el caso de una señal de alta frecuencia, la interferencia se produzca en espesores más pequeños. El siguiente fragmento de código crea una señal con frecuencias de 35 Hz, 55 Hz y 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]Un conjunto de señales fuente con frecuencias de 35 Hz, 55 Hz, 85 Hz.
Al calcular las trazas sísmicas y trazar gráficos de las amplitudes de las ondas reflejadas, podemos ver que para diferentes frecuencias se observa una anomalía en diferentes espesores de formación.
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]Gráficos de las amplitudes de la onda reflejada desde el borde superior de la cuña para diferentes frecuencias.
La capacidad de sacar conclusiones sobre el espesor de la formación a partir de los resultados de las observaciones sísmicas es extremadamente útil, porque una de las principales tareas en la exploración petrolera es evaluar los puntos más prometedores para colocar un pozo (es decir, aquellas áreas donde la formación está más grueso). Además, en la sección geológica puede haber objetos cuya génesis provoque un cambio brusco en el espesor de la formación. Esto hace que el análisis espectral sea una herramienta eficaz para estudiarlos. En la siguiente parte del artículo consideraremos estos objetos geológicos con más detalle.
Datos experimentales. ¿De dónde los conseguiste y qué buscar en ellos?
Los materiales analizados en el artículo se obtuvieron en Siberia occidental. La región, como probablemente todos, sin excepción, saben, es la principal región productora de petróleo de nuestro país. El desarrollo activo de los yacimientos comenzó en la región en los años 60 del siglo pasado. El principal método de búsqueda de depósitos de petróleo es la exploración sísmica. Es interesante observar imágenes de satélite de este territorio. A pequeña escala se puede observar una gran cantidad de pantanos y lagos; ampliando el mapa se pueden ver los sitios de perforación de pozos en racimo, y ampliando el mapa hasta el límite también se pueden distinguir los claros de los perfiles a lo largo de los cuales se producen fenómenos sísmicos. Se llevaron a cabo observaciones.
Imagen satelital de mapas Yandex - área de la ciudad de Noyabrsk
Una red de plataformas de pozos en uno de los campos.
Las rocas petrolíferas de Siberia occidental se encuentran en una amplia gama de profundidades, desde 1 km hasta 5 km. El mayor volumen de rocas que contienen petróleo se formó en el Jurásico y el Cretácico. El período Jurásico probablemente sea conocido por muchos por la película del mismo nombre. era significativamente diferente del moderno. La Enciclopedia Británica tiene una serie de paleomapas que caracterizan cada era helógica.
presente
Periodo Jurasico
Tenga en cuenta que en la época del Jurásico, el territorio de Siberia occidental era una costa marítima (tierra atravesada por ríos y un mar poco profundo). Dado que el clima era agradable, podemos suponer que un paisaje típico de aquella época era el siguiente:
Siberia jurásica
En esta imagen lo importante para nosotros no son tanto los animales y los pájaros, sino la imagen del río al fondo. El río es el mismo objeto geológico en el que nos detuvimos antes. El hecho es que la actividad de los ríos permite que se acumulen areniscas bien seleccionadas, que luego se convertirán en un depósito de petróleo. Estos embalses pueden tener una forma extraña y compleja (como el lecho de un río) y tienen un espesor variable: cerca de las orillas el espesor es pequeño, pero más cerca del centro del canal o en las zonas de meandros aumenta. Así, los ríos que se formaron en el Jurásico se encuentran ahora a una profundidad de unos tres kilómetros y son objeto de búsqueda de yacimientos de petróleo.
Datos experimentales. Procesamiento y visualización
Inmediatamente hagamos una reserva con respecto a los materiales sísmicos mostrados en el artículo, debido a que la cantidad de datos utilizados para el análisis es significativa, en el texto del artículo solo se incluye un fragmento del conjunto original de rastros sísmicos. Esto permitirá que cualquiera pueda reproducir los cálculos anteriores.
Cuando trabaja con datos sísmicos, un geofísico suele utilizar software especializado (hay varios líderes de la industria cuyos desarrollos se utilizan activamente, por ejemplo Petrel o Paradigm), que le permite analizar diferentes tipos de datos y tiene una interfaz gráfica conveniente. A pesar de toda la conveniencia, este tipo de software también tiene sus inconvenientes; por ejemplo, la implementación de algoritmos modernos en versiones estables lleva mucho tiempo y las posibilidades de automatizar los cálculos suelen ser limitadas. En tal situación, resulta muy conveniente utilizar sistemas matemáticos informáticos y lenguajes de programación de alto nivel, que permiten el uso de una amplia base algorítmica y, al mismo tiempo, asumen mucha rutina. Este es el principio utilizado para trabajar con datos sísmicos en Wolfram Mathematica. No es práctico escribir funciones ricas para el trabajo interactivo con datos; es más importante garantizar la carga desde un formato generalmente aceptado, aplicarles los algoritmos deseados y volver a cargarlos en un formato externo.
Siguiendo el esquema propuesto, cargaremos los datos sísmicos originales y los mostraremos en :
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]Los datos descargados e importados de esta forma son las rutas registradas en un área de 10 por 5 kilómetros. Si los datos se obtienen mediante un método de estudio sísmico tridimensional (las ondas no se registran a lo largo de perfiles geofísicos individuales, sino en toda el área simultáneamente), es posible obtener datos sísmicos cúbicos. Se trata de objetos tridimensionales, cuyas secciones verticales y horizontales permiten un estudio detallado del entorno geológico. En el ejemplo considerado, se trata de datos tridimensionales. Podemos obtener información del encabezado de texto, como este
StringPartition[seismic3DSEGY["textheader"], 80] // TableFormC 1 ESTE ES EL ARCHIVO DE DEMOSTRACIÓN PARA LA PRUEBA DEL PAQUETE GEOLOGYIO
C 2
C 3
C 4
C 5 FECHA NOMBRE DE USUARIO: USUARIO DE WOLFRAM
C 6 NOMBRE DE LA ENCUESTA: EN ALGÚN LUGAR DE SIBERIA
C 7 TIPO DE ARCHIVO VOLUMEN SÍSMICO 3D
C 8
C 9
GAMA C10 Z: PRIMEROS 2200M ÚLTIMOS 2400M
Este conjunto de datos será suficiente para que podamos demostrar las principales etapas del análisis de datos. Las huellas en el archivo se registran secuencialmente y cada una de ellas se parece a la siguiente figura: esta es la distribución de las amplitudes de las ondas reflejadas a lo largo del eje vertical (eje de profundidad).
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]Una de las trazas de la sección sísmica.
Sabiendo cuántas trazas hay en cada dirección del área estudiada, puede generar una matriz de datos tridimensional y mostrarla usando la función Image3D[]
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]]Imagen XNUMXD de un cubo de datos sísmicos (Eje vertical - profundidad)
Si las características geológicas de interés crean anomalías sísmicas intensas, entonces se pueden utilizar herramientas de visualización con transparencia. Las áreas "sin importancia" de la grabación se pueden hacer invisibles, dejando visibles sólo las anomalías. En Wolfram Mathematica esto se puede hacer usando и .
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]Imagen de cubo de datos sísmicos usando las funciones Opacity[] y Raster3D[]
Como en el ejemplo sintético, en secciones del cubo original se pueden identificar algunos límites geológicos (capas) con relieve variable.
La principal herramienta para el análisis espectral es la transformada de Fourier. Con su ayuda, puede evaluar el espectro de amplitud-frecuencia de cada traza o grupo de trazas. Sin embargo, después de transferir los datos al dominio de la frecuencia, se pierde información sobre en qué momentos (léase a qué profundidad) cambia la frecuencia. Para poder localizar los cambios de señal en el eje temporal (profundidad), se utilizan la transformada de Fourier en ventana y la descomposición wavelet. Este artículo utiliza la descomposición wavelet. La tecnología de análisis Wavelet comenzó a utilizarse activamente en la exploración sísmica en los años 90. Se considera que la ventaja sobre la transformada de Fourier en ventana es una mejor resolución temporal.
Usando el siguiente fragmento de código, puede descomponer una de las huellas sísmicas en componentes individuales:
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]Descomposición de una traza en componentes.
Para evaluar cómo se distribuye la energía de reflexión en diferentes tiempos de llegada de las ondas, se utilizan escalogramas (análogos a un espectrograma). Como regla general, en la práctica no es necesario analizar todos los componentes. Normalmente, se seleccionan componentes de frecuencia baja, media y alta.
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]Escalograma. Resultado de la función
Wolfram Language utiliza la función para la transformación wavelet . Y la aplicación de esta función a todo el conjunto de trazas se realizará mediante la función . Aquí vale la pena señalar uno de los puntos fuertes de Wolfram Mathematica: la capacidad de utilizar la paralelización. . En el ejemplo anterior, no hay necesidad de paralelización: el volumen de datos no es grande, pero cuando se trabaja con conjuntos de datos experimentales que contienen cientos de miles de rastros, esto es una necesidad.
tracesCWD=Table[Map[Hilbert[#,0]&,Re[ContinuousWaveletTransform[traces[[i]]][[1]]][[{13,15,18}]]],{i,1,Length@traces}]; Después de aplicar la función Aparecen nuevos conjuntos de datos correspondientes a las frecuencias seleccionadas. En el ejemplo anterior, estas frecuencias son: 38 Hz, 33 Hz, 27 Hz. La elección de las frecuencias se realiza con mayor frecuencia mediante pruebas: obtienen mapas efectivos para diferentes combinaciones de frecuencias y seleccionan la más informativa desde el punto de vista de un geólogo.
Si necesita compartir los resultados con colegas o proporcionárselos al cliente, puede utilizar la función SEGYExport[] del paquete 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];Con tres de estos cubos (componentes de baja frecuencia, media frecuencia y alta frecuencia), la combinación RGB se suele utilizar para visualizar los datos juntos. A cada componente se le asigna su propio color: rojo, verde, azul. En Wolfram Mathematica esto se puede hacer usando la función .
El resultado son imágenes a partir de las cuales se puede realizar una interpretación geológica. Los meandros que se registran en el tramo permiten delimitar paleocanales, que tienen más probabilidades de ser reservorios y contener reservas de petróleo. La búsqueda y el análisis de análogos modernos de dicho sistema fluvial nos permite determinar las partes más prometedoras de los meandros. Los propios canales se caracterizan por tener gruesas capas de arenisca bien seleccionada y son un buen depósito de petróleo. Las áreas fuera de las anomalías de "encaje" son similares a los depósitos de llanuras aluviales modernas. Los depósitos de las llanuras aluviales están representados principalmente por rocas arcillosas y la perforación en estas zonas será ineficaz.
Porción RGB del cubo de datos. En el centro (un poco a la izquierda del centro) se puede seguir el río serpenteante.
Porción RGB del cubo de datos. En el lado izquierdo se puede seguir el río serpenteante.
En algunos casos, la calidad de los datos sísmicos permite obtener imágenes mucho más claras. Esto depende de la metodología de trabajo de campo, del equipo utilizado por el algoritmo de reducción de ruido. En tales casos, no sólo son visibles fragmentos de sistemas fluviales, sino también paleoríos completos y extendidos.
Mezcla RGB de tres componentes de un cubo de datos sísmicos (corte horizontal). Profundidad aproximadamente 2 km.
Imagen de satélite del río Volga cerca de Saratov
Conclusión
Wolfram Mathematica le permite analizar datos sísmicos y resolver problemas aplicados relacionados con la exploración minera, y el paquete GeologyIO hace que este proceso sea más conveniente. La estructura de los datos sísmicos es tal que el uso de métodos integrados para acelerar los cálculos (, ,…) es muy eficiente y permite procesar grandes cantidades de datos. En gran medida, esto se ve facilitado por las funciones de almacenamiento de datos del paquete GeologyIO. Por cierto, el paquete se puede utilizar no sólo en el campo de la exploración sísmica aplicada. Casi los mismos tipos de datos se utilizan en radar de penetración terrestre y sismología. Si tiene sugerencias sobre cómo mejorar el resultado, qué algoritmos de análisis de señales del arsenal de Wolfram Mathematica son aplicables a dichos datos, o si tiene algún comentario crítico, por favor Deja un comentario.
Fuente: habr.com