Gracias al autor del blog Anton Ekimenko por su charla.Introduccion
Esta nota fue escrita a raíz de la
Conferencia de Tecnología Rusa de Wolfram y contiene una sinopsis del informe que di. El evento tuvo lugar en junio en la ciudad de San Petersburgo. Dado que trabajo a una cuadra del lugar de la conferencia, no pude evitar asistir a este evento. En 2016 y 2017, escuché los informes de la conferencia, y este año hice una presentación. En primer lugar, ha aparecido un tema interesante (como me parece a mí) que estamos desarrollando con
Kirill Belov , y en segundo lugar, después de un largo estudio de la legislación de la Federación de Rusia con respecto a la política de sanciones, han aparecido dos licencias de
Wolfram Mathematica en la empresa donde trabajo.
Antes de pasar al tema de mi presentación, me gustaría señalar la buena organización del evento. La página de la conferencia usa la imagen de la Catedral de Kazan. La Catedral es una de las principales atracciones de San Petersburgo y es muy claramente visible desde la sala en la que se celebró la conferencia.
A la entrada de la Universidad Estatal de Economía de San Petersburgo, los asistentes fueron recibidos por estudiantes asistentes, que no les permitieron perderse. Durante el registro, se entregaron pequeños recuerdos (un juguete con adhesivos, un bolígrafo, pegatinas con símbolos Wolfram). El almuerzo y las pausas para el café también se incluyeron en el horario de la conferencia. Sobre el delicioso café y los pasteles, ya noté en la pared del grupo: chefs bien hechos. Con esta parte introductoria, me gustaría enfatizar que el evento en sí, su formato y lugar ya traen emociones positivas.
El informe, que preparamos Kirill Belov y yo, se llama “Uso de Wolfram Mathematica para resolver problemas de geofísica aplicada. Análisis espectral de datos sísmicos o "dónde corrían los antiguos ríos". El contenido del informe abarca dos partes: en primer lugar, el uso de algoritmos disponibles en
Wolfram Mathematica para el análisis de datos geofísicos, y en segundo lugar, así es cómo colocar los datos geofísicos en Wolfram Mathematica.
Exploración sísmica
Primero debes hacer una pequeña excursión a la geofísica. La geofísica es una ciencia que estudia las propiedades físicas de las rocas. Bueno, dado que las rocas tienen diferentes propiedades: eléctricas, magnéticas, elásticas, existen métodos correspondientes de geofísica: exploración eléctrica, exploración magnética, exploración sísmica ... En el contexto de este artículo, solo tocaremos la exploración sísmica con más detalle. La exploración sísmica es el método principal para encontrar 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. La excitación de las vibraciones se realiza en tierra (dinamita o fuentes de vibración no explosivas de vibraciones elásticas) o en el mar (pistolas de aire comprimido). Las vibraciones elásticas se propagan a través del grosor de las rocas, siendo refractadas y reflejadas en los límites de las capas con diferentes propiedades. Las ondas reflejadas regresan a la superficie y son registradas por geófonos terrestres (típicamente dispositivos electrodinámicos basados en el movimiento de un imán suspendido en una bobina) o hidrófonos en el mar (según el efecto piezoeléctrico). En el momento de la llegada de las olas, uno puede juzgar las profundidades de las capas geológicas.
Equipo de remolque sísmico:
La pistola de aire excita las vibraciones elásticas:
Las olas pasan a través del macizo rocoso y son grabadas por hidrófonos:
Buque de investigación de exploración geofísica "Ivan Gubkin" en el muelle 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 las más importantes: la velocidad de propagación de las vibraciones elásticas y la densidad. Si dos capas tienen las mismas propiedades o propiedades similares, la onda no "notará" el límite entre ellas. Si las velocidades de onda en las capas difieren, entonces la reflexión ocurrirá en el límite de las capas. Cuanto mayor es la diferencia en las propiedades, más intensa es la reflexión. Su intensidad estará determinada por el coeficiente de reflexión (rc):
donde ρ es la densidad de la roca, ν es la velocidad de la onda, 1 y 2 denotan 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 se presenta una traza sísmica registrada como resultado de la convolución de una secuencia de coeficientes de reflexión con un pulso de sonda:
donde s (t) es la pista sísmica, es decir todo lo que un hidrófono o un geófono grabó durante un tiempo de grabación fijo,
w (t) es la señal que genera la pistola de aire,
n (t) es ruido aleatorio.
Calculamos, por ejemplo, una encuesta sísmica sintética. Utilizaremos el pulso Ricker ampliamente utilizado en la exploración sísmica como señal inicial.
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]
Fuente de pulso sísmico
Establecemos 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]]
La secuencia de coeficientes de reflexión:
Calcule y muestre la pista sísmica. Dado que los coeficientes de reflexión tienen signos diferentes, obtenemos dos reflexiones alternas en el camino sísmico.
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, se debe hacer una reserva: en realidad, la profundidad de las capas se determina, por supuesto, en metros, y el cálculo de la ruta sísmica se produce para el dominio del tiempo. Sería más correcto establecer las profundidades en metros y calcular los tiempos de llegada conociendo la velocidad en los estratos. En este caso, inmediatamente configuro las capas en el eje de tiempo.Si hablamos de estudios de campo, como resultado de tales observaciones, se registra una gran cantidad de tales series de tiempo (pistas sísmicas). Por ejemplo, al explorar un sitio de 25 km de largo y 15 km de ancho, donde, como resultado del trabajo, cada pista caracteriza una celda que mide 25x25 metros (dicha celda se llama bin), la matriz de datos final contendrá 600,000 pistas. Con un paso de muestreo de tiempo de 1 ms, un tiempo de grabación de 5 segundos, el archivo de datos final tendrá más de 11 GB, y la cantidad de materia prima original puede ser de cientos de gigabytes.
¿Cómo trabajar con ellos en
Wolfram Mathematica ?
El desarrollo del paquete comenzó con una
pregunta 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. La
publicación correspondiente
en el muro de Wolfram Community incluso
fue notada por los moderadores. Actualmente, el paquete admite trabajar con los siguientes tipos de datos que se utilizan activamente en la industria geológica:
- importar datos de mapas en formato ZMAP e IRAP
- importación de mediciones en pozos de formato LAS
- entrada y salida de archivos sísmicos en formato SEGY
Para instalar el paquete, debe seguir las instrucciones en la página de descarga del paquete ensamblado, es decir, ejecute el siguiente código en cualquier
cuaderno de Mathematica :
If[PacletInformation["GeologyIO"] === {}, PacletInstall[URLDownload[ "https://wolfr.am/FiQ5oFih", FileNameJoin[{CreateDirectory[], "GeologyIO-0.2.2.paclet"}] ]]]
Después de eso, el paquete se instalará en la carpeta predeterminada, cuya ruta se puede obtener de la siguiente manera:
FileNameJoin[{$UserBasePacletsDirectory, "Repository"}]
Como ejemplo, demostramos las características principales del paquete. La convocatoria es tradicionalmente para paquetes en Wolfram Language:
Get["GeologyIO`"]
El paquete se está desarrollando con
Wolfram Workbench . Esto le permite acompañar la funcionalidad principal del paquete con la documentación, que no difiere en formato de presentación de la documentación de Wolfram Mathematica y proporciona al paquete los archivos de prueba para el primer conocido.
Tal archivo, en particular, es el archivo Marmousi.segy, un modelo sintético de una sección geológica desarrollado por el Instituto Francés del Petróleo. Con este modelo, los desarrolladores prueban sus propios algoritmos para modelar el campo de onda, el procesamiento de datos, la inversión de trazas sísmicas, etc. El modelo Marmousi mismo 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"]
El resultado de la importación es un objeto SEGYData:
El formato SEGY implica el almacenamiento de diversa información sobre las observaciones. En primer lugar, estos son comentarios de texto. Aquí se ingresa información sobre el lugar de trabajo, 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 Marmousi se generó en el Instituto ... velocidad mínima de 1500 m / sy un máximo de 5500 m / s)"
El modelo geológico real puede mostrarse contactando las huellas sísmicas utilizando la tecla "trazas" (una de las características del paquete es la insensibilidad a mayúsculas y minúsculas de las teclas):
ArrayPlot[Transpose[marmousi["traces"]], PlotTheme -> "Detailed"]
Modelo Marmousi:
Por el momento, el paquete también le permite descargar datos en partes de archivos grandes, para que sea posible procesar archivos cuyo tamaño puede alcanzar decenas de gigabytes. El paquete también incluye funciones para exportar datos a .segy y reescritura parcial al final del archivo.
Por separado, vale la pena señalar la funcionalidad del paquete cuando se trabaja con la compleja estructura de los archivos .segy. Dado que permite no solo acceder a claves e índices para rastreos individuales, encabezados, sino también cambiarlos con posterior escritura en un archivo. Muchos de los detalles técnicos de la implementación de GeologyIO están más allá del alcance de este artículo y probablemente merecen una descripción por separado.
La relevancia del análisis espectral en la exploración sísmica.
La capacidad de importar materiales sísmicos en Wolfram Mathematica le permite utilizar la funcionalidad de procesamiento de señal incorporada para datos experimentales. Como cada pista sísmica es una serie temporal, una de las principales herramientas para estudiarlas es el análisis espectral. Entre los requisitos previos para el análisis de la composición de frecuencia de los datos sísmicos se encuentran, por ejemplo, los siguientes:
- Los diferentes tipos de ondas se caracterizan por una composición de frecuencia diferente. Esto nos permite seleccionar ondas útiles y suprimir las ondas de interferencia.
- Las propiedades de la roca, como la porosidad y la saturación, pueden afectar la composición de frecuencia. Esto le permite seleccionar rocas con mejores propiedades.
- Las capas con diferentes espesores causan anomalías en diferentes rangos de frecuencia.
El tercer párrafo es central en el contexto de este artículo. A continuación se muestra un fragmento de código para calcular trazas sísmicas en el caso de una capa con grosor variable: un modelo de cuña. Este modelo se estudia tradicionalmente en la exploración sísmica para el análisis de los efectos de interferencia, cuando las ondas reflejadas desde 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 de arrebato:
La velocidad de la onda dentro de la cuña es de 4500 m / s, fuera de la cuña es de 4000 m / s, y se supone que la densidad es constante de 2200 g / cm³. Para tal modelo, calculamos los coeficientes de reflexión y las huellas 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]
Pistas sísmicas para el modelo de cuña:
La secuencia de huellas sísmicas que se muestra en esta figura se denomina sección sísmica. Como puede ver, su interpretación se puede realizar en un nivel intuitivo, ya que la geometría de las ondas reflejadas corresponde únicamente al modelo que se estableció anteriormente. Si analiza las pistas con más detalle, puede ver que las pistas del 1 al 30 no difieren: el reflejo desde la parte superior de la formación y desde la suela no se superponen. A partir de la ruta 31, los reflejos comienzan a interferir. Y, aunque, en el modelo, los coeficientes de reflexión no cambian horizontalmente, las pistas sísmicas cambian su intensidad con un cambio en el grosor del depósito.
Considere la amplitud de la reflexión desde el límite superior del depósito. A partir de la ruta 60, la intensidad de la reflexión comienza a aumentar y en la ruta 70 se vuelve máxima. Así es como se manifiesta la interferencia de las olas desde el techo y la suela para los estratos, lo que lleva en algunos casos a anomalías significativas 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 frecuencia más baja, entonces la interferencia comienza a aparecer en grandes espesores de la formación, y en el caso de una señal de alta frecuencia, la interferencia ocurre 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
Después de calcular los rastros sísmicos y trazar las amplitudes de la onda reflejada, podemos ver que para diferentes frecuencias se observa una anomalía en diferentes espesores de la 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 grosor del reservorio en función de los resultados de las observaciones sísmicas es extremadamente útil, porque una de las principales tareas en la prospección de un campo petrolero es evaluar los puntos más prometedores para la colocación de un pozo (es decir, aquellas áreas donde el reservorio tiene un gran grosor). Además, en la sección geológica, puede haber objetos que, debido a su génesis, causen un cambio brusco en el grosor de la formación. Esto hace que el análisis espectral sea una herramienta efectiva para estudiarlos. En la siguiente parte del artículo, consideraremos tales objetos geológicos con más detalle.
Datos experimentales ¿Dónde se reciben y qué buscar en ellos?
Los materiales que se analizan en el artículo se obtuvieron en el territorio de Siberia occidental. La región, como probablemente todos lo saben, es la principal región productora de petróleo de nuestro país. El desarrollo activo de los nacimientos en el metro comenzó en la región en los años 60 del siglo pasado. El método principal para encontrar campos petroleros es la exploración sísmica. Es interesante considerar imágenes satelitales de este territorio. A pequeña escala, se puede observar una gran cantidad de pantanos y lagos, al aumentar el mapa, puede ver los pozos de conglomerados para perforar pozos, y al aumentar el mapa hasta el límite, también puede distinguir los claros de los perfiles mediante los cuales se realizaron observaciones sísmicas.
Imagen satelital de los mapas de Yandex - el área de la ciudad de Noyabrsk:
Una red de sitios de clúster en uno de los campos:
Las rocas petroleras de la Siberia occidental se encuentran en una amplia gama de profundidades, desde 1 km hasta 5 km. El volumen principal de rocas que contienen petróleo se forma en los tiempos jurásico y cretáceo. El período jurásico es probablemente conocido por muchos por la película del mismo nombre.
El clima del período Jurásico fue significativamente diferente del moderno. En la Enciclopedia Británica hay una serie de paleocartas que caracterizan cada era geológica.
Tiempo presente:
Periodo jurásico:
Tenga en cuenta que en el Jurásico, el territorio de Siberia occidental era una costa marítima (tierra atravesada por ríos y mar poco profundo). Como el clima era cómodo, se puede suponer que un paisaje típico de esa época tenía el siguiente aspecto:
Siberia del Jurásico:
En esta imagen, no tanto los animales y las aves son importantes para nosotros, como la imagen del río en el fondo. Un río es el objeto muy geológico en el que nos detuvimos antes. El hecho es que la actividad de los ríos permite la acumulación de areniscas bien clasificadas, que luego se convierten en un depósito de petróleo. Estos depósitos pueden tener una forma extraña y compleja (como el lecho de un río) y tienen un grosor variable: a lo largo de la costa, el grosor es pequeño y aumenta más cerca del centro del canal o en secciones de los meandros. Entonces, los ríos formados en el tiempo jurásico ahora están a una profundidad de aproximadamente tres kilómetros y son el objeto de la búsqueda de depósitos de petróleo.
Datos experimentales Procesamiento y visualización.
Inmediatamente hacemos una reserva con respecto a los materiales sísmicos que se muestran en el artículo, debido al hecho de que la cantidad de datos utilizados para el análisis es significativa, solo un fragmento del conjunto original de rastros sísmicos se coloca en el texto del artículo. Esto permitirá a todos reproducir los cálculos anteriores.
Cuando se trabaja con datos sísmicos, los geofísicos generalmente usan software especializado (hay varios líderes de la industria cuyo desarrollo se usa 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 comodidad, 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 automatización de 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 le permitan utilizar una amplia base algorítmica y, al mismo tiempo, asumir mucha rutina. Utilizando este principio, se construye el trabajo con datos sísmicos en Wolfram Mathematica. No es práctico escribir una gran cantidad de trabajo interactivo con datos; es más importante garantizar la carga desde un formato generalmente aceptado, aplicarles los algoritmos deseados y volver a cargarlo en un formato externo.
Siguiendo el esquema propuesto, cargue los datos sísmicos originales y muéstrelos en
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]
Los datos descargados e importados de esta manera son las pistas registradas en un sitio que mide 10 por 5 kilómetros. En el caso de que los datos se obtengan por el método de exploración sísmica tridimensional (las ondas no se registran a lo largo de perfiles geofísicos separados, sino en toda el área al mismo tiempo), es posible obtener cubos de datos sísmicos. Estos son objetos tridimensionales, cuyas secciones verticales y horizontales permiten un estudio detallado del entorno geológico. En el ejemplo considerado, estamos tratando solo con datos tridimensionales. Podemos obtener información del encabezado del texto, por ejemplo, así
StringPartition[seismic3DSEGY["textheader"], 80] // TableForm
C 1 ESTE ES EL ARCHIVO DEMO PARA LA PRUEBA DE PAQUETE GEOLOGYIO
C 2
C 3
C 4
C 5 FECHA NOMBRE DEL USUARIO: USUARIO DE WOLFRAM
NOMBRE DE LA ENCUESTA C 6: EN ALGÚN LUGAR DE SIBERIA
C 7 TIPO DE ARCHIVO 3D VOLUMEN SÍSMICO
C 8
C 9
GAMA C10 Z: PRIMEROS 2200M ÚLTIMOS 2400M
Este conjunto de datos será suficiente para que demostremos las etapas principales del análisis de datos. Las pistas en el archivo se graban 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 secciones sísmicas:
Sabiendo cuántos rastros se encuentran en cada dirección del área estudiada, puede formar una matriz de datos tridimensionales 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 tridimensional de un cubo de datos sísmicos. (Eje vertical - profundidad):
En el caso de que los objetos geológicos de interés creen anomalías sísmicas intensas, se pueden utilizar herramientas de visualización con transparencia. Las secciones "sin importancia" de la grabación pueden hacerse invisibles, dejando solo anomalías visibles. En Wolfram Mathematica, esto se puede hacer usando
Opacity [] y
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]
Imagen de un cubo de datos sísmicos utilizando las funciones Opacidad [] y Raster3D []
Como en el ejemplo sintético, en las rebanadas del cubo original, se pueden distinguir algunos límites geológicos (capas) con relieve variable.
La principal herramienta de análisis espectral es la transformada de Fourier. Utilizándolo, puede evaluar el espectro de amplitud-frecuencia de cada ruta o grupo de rutas. Sin embargo, después de transferir datos al dominio de frecuencia, se pierde información sobre a qué horas (leer a qué profundidades) cambia la frecuencia. Para poder localizar los cambios de señal en el eje de tiempo (profundo), se utilizan la transformación de Fourier de ventana y la descomposición de wavelet. Este artículo utiliza la descomposición wavelet. La tecnología de análisis Wavelet comenzó a usarse activamente en la exploración sísmica en los años 90. La ventaja sobre la ventana de la transformada de Fourier se considera la mejor resolución de tiempo.
Usando el siguiente fragmento de código, puede descomponerse en componentes individuales de uno de los rastros sísmicos:
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 componentes
Para evaluar cómo se distribuye la energía de reflexión en diferentes momentos de llegada de las ondas, se utilizan escalogramas (análogos del espectrograma). Como regla general, en la práctica, no es necesario analizar todos los componentes. Por lo general, elija el componente de baja, media y alta frecuencia.
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. Resultado de la función
WaveletScalogram []
Wolfram Language utiliza la función ContinuousWaveletTransform [] para las transformaciones wavelet . Y la aplicación de esta función a todo el conjunto de trazas se lleva a cabo con el uso de la función Tabla [] . Aquí vale la pena señalar que uno de los puntos fuertes de Wolfram Mathematica es la capacidad de usar la paralelización de ParallelTable [] . En el ejemplo dado, la paralelización no es necesaria: el volumen de datos no es grande, pero cuando se trabaja con conjuntos de datos experimentales que contienen cientos de miles de trazas, 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 ContinuousWaveletTransform [] , aparecen nuevas matrices de datos correspondientes a las frecuencias seleccionadas. En el ejemplo anterior, estas son frecuencias: 38Hz, 33Hz, 27Hz. La elección de frecuencias se realiza con mayor frecuencia en base a pruebas: obtienen mapas efectivos para diferentes combinaciones de frecuencias y eligen 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 usar la función SEGYExport [] de 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];
Teniendo tres de estos cubos disponibles (componentes de baja frecuencia, frecuencia media y alta frecuencia), por regla general, la mezcla RGB se usa para la visualización conjunta de datos. A cada componente se le asigna su propio color: rojo, verde, azul. En Wolfram Mathematica, esto se puede hacer usando la función ColorCombine [] .Como resultado, se obtienen imágenes que pueden usarse para realizar una interpretación geológica. Los meandros que cuelgan de la rebanada le permiten contornear el paleorus, que es más probable que sean reservorios y contengan 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 lechos en sí se caracterizan por gruesas capas de arenisca bien clasificada y son un buen depósito para el petróleo. Las áreas fuera de las anomalías de "encaje" son similares a los depósitos modernos de llanuras aluviales. Los sedimentos de la llanura de inundación están representados principalmente por rocas arcillosas y la perforación en estas zonas será ineficaz.Rebanada RGB del cubo de datos. En el centro (algo a la izquierda del centro) puedes rastrear el río serpenteante.
Rebanada RGB del cubo de datos. En el lado izquierdo se puede rastrear el río serpenteante.
En algunos casos, la calidad de los datos sísmicos permite imágenes significativamente más nítidas. Depende de la metodología de trabajo de campo, equipo adoptado por el algoritmo de reducción de ruido. En tales casos, no solo son visibles fragmentos del sistema fluvial, sino también paleorecs extendidos enteros.Mezcla RGB de los tres componentes del cubo de datos sísmicos (corte horizontal). La profundidad es de aproximadamente 2 km.
Imagen satelital del río Volga en la región de Saratov
Conclusión
En Wolfram Mathematica, puede analizar datos sísmicos y resolver problemas de aplicación asociados con la búsqueda de minerales, y el paquete GeologyIO le permite hacer este proceso más conveniente. La estructura de los datos sísmicos es tal que el uso de métodos integrados para acelerar los cálculos ( ParallelTable [] , ParallelDo [], ...) es muy eficiente y le permite procesar grandes cantidades de datos. En gran medida, esto se ve facilitado por las características de almacenamiento de datos del paquete GeologyIO. Por cierto, el paquete puede usarse no solo en el campo de la exploración sísmica aplicada. Se utilizan casi los mismos tipos de datos en georadar y sismología. Si tiene sugerencias sobre cómo mejorar el resultado, qué algoritmos de análisis de señal de Wolfram Mathematica son aplicables a dichos datos, o si tiene comentarios críticos, deje comentarios.