A proposal for regularized inversion for an ill-conditioned deconvolution operator

Desde el punto de vista de la teoría de problemas inversos, la deconvolución puede ser entendida como una inversión lineal de un problema mal-puesto y mal-condicionado. La característica del mal-condicionamiento del operador de deconvolución hace que la solución del problema inverso sea sensitiva a errores en los datos. La regularización de Tikhonov es el método más común empleado para estabilizar la solución y obtener su unicidad. Sin embargo, los resultados del método de Tikhonov no proveen calidad suficiente cuando el ruido en los datos es fuerte. Este trabajo hace uso del método del gradiente conjugado, basado en el esquema de Tikhonov aplicado a la deconvolución, cuyo parámetro de regularización es calculado iterativamente teniendo en cuenta el criterio de discrepancia de Morozov en la función objetivo. Haciendo uso de datos sísmicos sintéticos como datos reales apilados, se llevó a cabo el proceso de deconvolución con y sin regularización basada en el algoritmo del gradiente conjugado. Se llevó a cabo una comparación del esquema planteado. Aplicando la deconvolución regularizada en los datos sintéticos muestra una mejora en la estabilidad de la solución y los datos sísmicos post-apilados mostraron un incremento de la resolución vertical aun con ruido en los datos.From the inverse problem theory aspect, deconvolution can be understood as the linear inversion of an ill-posed and ill-conditioned problem. The ill-conditioned property of the deconvolution operator make the solution of inverse problem sensitive to errors in the data. Tikhonov regularization is the most commonly used method for stability and uniqueness of the solution. However, results from Tikhonov method do not provide sufficient quality when the noise in the data is strong. This work uses the conjugate gradient method applied to the Tikhonov deconvolution scheme, including a regularization parameter calculated iteratively and based on the improvement criterion of Morozov discrepancy applied on the objective function. Using seismic synthetic data and real stacked seismic data, we carried out a deconvolution process with regularization and without regularization based on a conjugated gradient algorithm. A comparison of results is also presented. Applying regularized deconvolution on synthetic data shows improved stability on the solution. Additionally, real post-stack seismic data shows a direct application for increasing the vertical resolution even with noisy data.Desde o ponto de vista da teoria de problemas inversos, a deconvolução pode ser entendida como uma inversão linear de um problema mal posto e mal condicionado. A característica do mal condicionamento do operador de deconvolução faz que a solução do problema inverso seja sensitiva a erros nos dados. A regularização de Tikhonov é o método mais comum utilizado para estabilizar a solução e obter sua unicidade. Porém, os resultados do método de Tikhonov não fornecem qualidade suficiente quando o ruído nos dados é alto. Este trabalho faz uso do método do gradiente conjugado, baseado no esquema de Tikhonov aplicado à deconvolução, cujo parâmetro de regularização é calculado iterativamente tendo em conta o critério de discrepância de Morozov na função objetivo. Fazendo uso de dados sísmicos sintéticos como dados reais empilhados, foi realizado o processo de deconvolução com e sem regularização baseado no algoritmo do gradiente conjugado. Realizou-se uma comparação do esquema proposto. Aplicando a deconvolução regularizada nos dados sintéticos mostra uma melhoria na estabilidade da solução e os dados sísmicos pós-empilhados mostraram um aumento da resolução vertical mesmo com ruído nos dados

Deconvolution type imaging condiction effects on short profile migration amplitudes

Amplitude preservation in Pre-Stack Depth Migration (PSDM) processes that use wave field extrapolation must be ensured first, in the operators used to continue the wavefield in time or depth, and second, in the imaging condition used to estimate the reflectivity function. In the later point, the conventional correlation-type imaging condition must be replaced by a deconvolution-type imaging condition. Migration performed in common-shot profile domain obtains the final migrated image as the superposition of images resulting of migrate each shot separately. The amplitude obtained in a point of the migrated image corresponds to the sum of the reflectivities for each shot which has illuminated such point, along the angles determined by the velocity model and the positions of the source and the receiver. The deeper the reflector, the lower the amplitude of the illumination field will be. As result, the correlation-type imaging condition produces images with an unbalanced amplitude decrease with depth. A deconvolution-type imaging condition scales the amplitudes through a correlation, using the weighting function dependent on the spectral density or the illumination of the downgoing wave field. In this article, two possible scaling functions have been used in the case of a single shot. In the case of data with multiple shots, five scaling possibilities are presented with the spectral density or the illumination function. The results of applying these imaging conditions to synthetic data with multiple shots show that the values of the amplitude in the migrated images are influenced by the coverage of the common midpoint, compensating this effect only in one of the imaging conditions described. Numerical experiments with synthetic data generated using Seismic Unix and the Sigsbee2a data are presented, highlighting that in velocity fields with strong vertical and lateral velocity variations, the balance of the amplitudes of the deep reflectors relative to the shallow reflectors is strongly influenced by the imaging condition applied.A preservação de amplitudes nos processos PreStack Depth Migration (PSDM pelas suas siglas em inglês) que utilizam extrapolação de campos de onda, deve-se garantir tantonos operadores que são utilizados para continuar os campos em tempo ou profundidade, como na condição de imagem que é utilizada para estimar a função de refletividade, para o qual a convencional condição de imagem de correlação deve ser substituída por uma condição de imagem tipo deconvolução. A migração realizada no domínio de disparo comum obtém a imagem migrada final, como a superposição das imagens obtidas de migrar os disparos separadamente. A amplidão obtida em um ponto da imagem migrada corresponde à soma das refletividades obtidas para cada disparo, que alcança a iluminar dito ponto a um ângulo determinado pelo modelo de velocidades e as posições da fonte e o receptor. Quanto mais profundos os refletores, menores serão as amplitudes do campo de iluminação e, portanto, uma condição de imagem tipo correlação sempre produzirá imagens com amplitudes desbalanceadas que diminuem com a profundidade. A condição de imagem tipo deconvolução busca realizar um escalamento das amplitudes obtidas através da correlação, utilizando una função de peso dependente da densidade espectral ou a iluminação do campo descendente. No caso de um único disparo, duas possibilidades de escalamento foram utilizadas. No caso de dados com múltiplos disparos se apresentam cinco possibilidades de escalamento com a densidade espectral ou a função iluminação. Os resultados da aplicação destas condições de imagem a dados sintéticos com múltiplos disparos mostram que os valores de amplidão nas imagens migradas estão influenciados pela cobertura do ponto médio comum, sendo compensado este efeito somente em uma das condições de imagem descritas. Experimentos numéricos com dados sintéticos gerados em Seismic Unix e com o dado Sigsbee2a são apresentados, ressaltando que em modelos com campos de velocidades que incluem fortes variações verticais e laterais de velocidade, o balanço das amplitudes dos refletores profundos em relação aos refletores rasos é fortemente dependente da condição de imagem aplicada.La preservación de amplitudes en los procesos PreStack Depth Migration (de las siglas en inglés PSDM) que utilizan extrapolación de campos de onda, se debe garantizar tanto en los operadores que son utilizados para continuar los campos en tiempo o profundidad, como en la condición de imagen que es utilizada para estimar la función de reflectividad, para lo cual la convencional condición de imagen de correlación debe ser remplazada por una condición de imagen tipo deconvolución. La migración realizada en el dominio de disparo común, obtiene la imagen migrada final como la superposición de las imágenes obtenidas de migrar los disparos separadamente. La amplitud obtenida en un punto de la imagen migrada corresponde a la suma de las reflectividades obtenidas para cada disparo que alcanza a iluminar dicho punto, a un ángulo determinado por el modelo de velocidades y las posiciones de la fuente y el receptor. Cuanto más profundos los reflectores, menores serán las amplitudes del campo de iluminación y por lo tanto una condición de imagen tipo correlación siempre producirá imágenes con amplitudes desbalanceadas que disminuyen con la profundidad. La condición de imagen tipo deconvolución busca realizar un escalamiento de las amplitudes obtenidas a través de la correlación, utilizando una función peso dependiente de la densidad espectral o la iluminación del campo descendente. En el caso de un único disparo, dos posibilidades de escalamiento han sido utilizadas. En el caso de datos con múltiples disparos se presentan cinco posibilidades de escalamiento con la densidad espectral o la función iluminación. Los resultados de la aplicación de estas condiciones de imagen a datos sintéticos con múltiples disparos muestran que los valores de amplitud en las imágenes migradas están influenciados por la cobertura del punto medio común, siendo compensado este efecto solamente en una de las condiciones de imagen descritas. Experimentos numéricos con datos sintéticos generados en Seismic Unix y con el dato Sigsbee2a son presentados, resaltando que en modelos con campos de velocidades que incluyen fuertes variaciones verticales y laterales de velocidad, el balance de las amplitudes de los reflectores profundos en relación a los reflectores rasos es fuertemente dependiente de la condición de imagen aplicada

A practical implementation of acoustic full waveform inversion on graphical processing units

Recently, Full Waveform Inversion (FWI) has gained more attention in the exploration geophysics community as a data fitting method that provides high-resolution seismic velocity models. Some of FWI essential components are a cost function to measure the misfit between observed and modeled data, a wave propagator to compute the modeled data and an initial velocity model that is iteratively updated until an acceptable decrease of the cost function is reached. Since FWI is a wave equation based method, the computational costs are elevated. In this paper, it is presented a fast Graphical Processing Unit (GPU) FWI implementation that uses a 2D acoustic wave propagator in time and updates the model using the gradient of the cost function, which is efficiently computed with the adjoint state method. The proposed parallel implementation is tested using the Marmousi velocity model. The performance of the proposed implementation is evaluated using the NVIDIA GeForce GTX 860 GPU and compared to a serial Central Processing Unit (CPU) implementation, in terms of execution time. We also evaluate the GPU occupancy and analyze the memory requirements. Our tests show that the GPU implementation can achieve a speed-up of 26.89 times when compared to its serial CPU implementation.Recentemente, a inversão de onda completa (FWI, sigla em inglês) ganhou maior atenção na comunidade de exploração geofísica como método de ajuste de dados, que fornece modelos de velocidades sísmicas de alta resolução. Alguns dos componentes essenciais do FWI são uma função de custo para estimar a diferença entre os dados observados e os dados modelados, um propagador do campo de ondas acústicas para os dados modelados e um modelo de velocidade inicial, que é atualizada de forma iterativa. Como o FWI está baseado no método da equação da onda, as exigências computacionais de execução são altas. Neste artigo apresentamos uma implementação rápida do FWI acústico 2D em tempo em uma unidade de processamento gráfico (GPU, sigla em inglês). Esta implementação utiliza um propagador da equação de onda e atualiza o modelo de velocidade, utilizando o gradiente da função objetivo, que é calculada de forma eficiente usando o método do estado adjunto. Proposta de implementação paralela é testada utilizando o modelo de velocidade Marmousi. O desempenho da implementação proposta é avaliada usando uma GeForce GTX 860 e comparada com uma aplicação de série em, um único processador, em termos de tempo de execução. Avaliamos também a quantidade de recursos utilizados pela GPU e analisamos os requisitos de memória. Os testes mostram que a implementação em GPU pode conseguir uma taxa de aceleração de 26.89 vezes quando comparada com uma implementação serial do processador.Recientemente, la inversión de onda completa (FWI, por sus siglas en inglés) ha ganado una mayor atención en la comunidad de exploración geofísica como un método de ajuste de datos que provee modelos de velocidades sísmicas de gran resolución. Algunos de los componentes esenciales de la FWI corresponden a una función de costo para medir la diferencia entre los datos observados y los datos modelados, un propagador de onda para obtener los datos modelados y un modelo de velocidad inicial que es actualizado iterativamente hasta llegar a un valor deseado de la función de costo. Como la FWI es un método basado en la ecuación de onda, el costo computacional de su implementación es elevado. En este documento presentamos una implementación rápida de la FWI 2D acústica en tiempo sobre una unidad de procesamiento gráfico (GPU, por sus siglas en inglés). Esta implementación usa la ecuación de onda acústica para modelar la propagación y actualiza el modelo de velocidades usando el gradiente de la función de costo, el cual es calculado eficientemente usando el Método del Estado Adjunto. La implementación paralela propuesta es probada usando el modelo de velocidades Marmousi. El desempeño de la implementación propuesta es evaluado usando una GPU NVIDIA GeForce GTX 860 y comparado con una implementación serial sobre un procesador, en términos de tiempo de ejecución. Adicionalmente, se evalúa la cantidad de recursos usados por la GPU y se analizan los requerimientos de memoria de la implementación. Las pruebas muestran que la implementación sobre GPU puede alcanzar un índice de aceleración de 26.89 veces si se compara con la implementación serial sobre el procesador