Stabilized inverse Q filtering

Stabilized inverse Q filtering is a data processing technology for enhancing the resolution of reflection seismology images where the stability of the method used is considered. Q is the anelastic attenuation factor or the seismic quality factor, a measure of the energy loss as the seismic wave moves. To obtain a solution when we make computations with a seismic model we always have to consider the problem of instability and try to obtain a stabilized solution for seismic inverse Q filtering.

When a wave propagates through subsurface materials both energy dissipation and velocity dispersion takes place. Inverse Q filtering is a method to restore the energy loss due to energy dissipation (amplitude compensation) and to correct the time-shift of the data due to velocity dispersion.

Wang has written an excellent book on the subject of inverse Q filtering, Seismic inverse Q filtering (2008), and discuss the subject of stabilizing the method. He writes: “The phase-only inverse Q filter mentioned above is unconditionally stable. However, if including the accompanying amplitude compensation in the inverse Q filter, stability is a major issue of concern in implementation.”^[1]

Hale (1981)^[2] found that the inverse Q filter overcompensated the amplitudes for the later events in a seismic trace. Therefore, in order to obtain reasonable amplitude, the amplitude spectrum of the computed filter has to be clipped at some maximum gain to prevent undue amplitude at later times. On basis of this concept Wang proposed a stabilized inverse Q filtering approach that was able to compensate simultaneously for both attenuation and dispersion.”^[3] The unclipped version of Wang’s solution is presented in the wikipedia article seismic inverse Q filtering. The solution is based on the theory of wavefield downward continuation. In this outline here I will compute on a clipped version by introducing low-pass filtering. Both Hale and Wang introduced low-passfiltering as a method for stabilization.

We have the equation for seismic inverse Q filtering from Wang:

${\displaystyle U(t+\bigtriangleup t,w)=U(t,w)\exp {\bigg (}|{\frac {w}{w_{r}}}|^{-\gamma }{\frac {|w|\bigtriangleup t}{2Q(w)}}{\bigg )}\exp {\bigg (}i|{\frac {w}{w_{r}}}|^{-\gamma }w\bigtriangleup t{\bigg )}\quad (1)}$

${\displaystyle \gamma =(\pi Q_{r})^{-1}}$

Time is denoted τ, frequency is w and i is the imaginary unit. Q_r and w_r are reference values representing damping and frequency for a certain frequency. To demonstrate stability we can simply bypass using a reference frequency and get a more simple equation:

${\displaystyle U(t+\bigtriangleup t,w)=U(t,w)\exp {\bigg (}{\frac {|w|\bigtriangleup t}{2Q(w)}}{\bigg )}\exp {\bigg (}iw\bigtriangleup t{\bigg )}\quad (2)}$

The sum of these plane waves gives the time-domain seismic signal,

${\displaystyle U(t+\bigtriangleup t)=\int _{0}^{\infty }U(t+\bigtriangleup t,w)dw.\quad (2.b)}$

En la figura 1 se presenta la solución de (2/2.b) para un modelo sísmico con diferentes valores de Q, lo que evidencia claramente la inestabilidad numérica. El número en la parte superior de la figura 1 corresponde al valor de Q: 1=Q1, 2=Q2, etc. Los resultados son similares a los presentados en el libro de Wang (cada traza se escala individualmente, por lo que los artefactos son más pronunciados en la traza 5 que en la traza 4). Sin embargo, Wang también consideró la compensación de fase. Los cálculos aquí realizados corresponden únicamente a la inversión de amplitud, ya que la compensación de fase no es necesaria para demostrar la inestabilidad, puesto que siempre es estable.

• 
  Figura 1. Trazas filtradas con Q inverso Q1=400,Q2=200,Q3=100,Q4=50,Q5=25

En la práctica, los artefactos causados ​​por la inestabilidad numérica pueden suprimirse mediante un filtro de paso bajo. Hale escribió que el IQF sin recorte de un sismograma amplificaba la frecuencia de Nyquist por un factor de 7 millones cuando teníamos la relación t/Q=10 y concluyó que para sismogramas típicos con longitudes mayores a 1000 muestras y un valor Q alrededor de 100, los datos rara vez son lo suficientemente puros como para justificar el uso del IQF sin recorte. Wang introdujo una frecuencia de corte para establecer un criterio para la estabilización mediante una fórmula matemática. Sin embargo, considerando el artículo de Hale, podría ser suficiente con simplemente eliminar la frecuencia de Nyquist. Esto significa que la frecuencia cercana a la frecuencia de Nyquist sea la frecuencia de corte. En la figura 2 vemos un modelo sísmico que nos proporciona datos de referencia para el filtrado Q inverso (gráfico rojo). Veremos que el IQF de este modelo amplificará la frecuencia de Nyquist por un factor ligeramente inferior a 5 millones.

• 
  Figura 2. Modelo de datos de referencia sísmica. El gráfico verde representa la respuesta al impulso sin amortiguación y el gráfico rojo la respuesta al impulso amortiguada Q=50.

La figura 3 muestra la traza filtrada mediante Q inversa, considerando solo la amplitud, de la figura 2 para Q=50 (traza 4). El resultado evidencia la inestabilidad numérica. Se observan artefactos a lo largo de toda la traza.

• 
  Fig. 3. Modelo de la fig. 2. con filtro Q inverso aplicado Q=50.

Intentaremos eliminar los artefactos aplicando un filtro de paso bajo a la traza de la figura 3. Utilizamos la herramienta de procesamiento de señales de MATLAB y creamos un filtro de paso bajo (filtro IIR de fase cero) en la figura 4 con una frecuencia de corte de 120  Hz. La respuesta de amplitud del filtro se muestra en azul y la fase en verde.

• 
  Figura 4. Filtro de paso bajo

El resultado de filtrar la traza de la figura 3 con el filtro de paso bajo de la figura 4 se muestra en la figura 5. Se eliminan todos los artefactos y nos queda la respuesta impulsional que se puede comparar con el modelo original de la figura 2.

• 
  Figura 5. Respuestas impulsionales filtradas donde se elimina toda la energía inestable.

Un estudio de la respuesta en frecuencia de la traza de la figura 3 (sin recorte) y la figura 5 (recortada) proporcionará una mejor comprensión del proceso de filtrado. La figura 6 muestra la magnitud de la respuesta en frecuencia en función de la frecuencia digital antes del filtrado. Esta representación ofrece una buena imagen de lo que sucede alrededor de la frecuencia de Nyquist cuando se aplica el filtro paso bajo. Se acumula energía inestable cerca de la frecuencia de Nyquist. Después del filtrado, la energía inestable alrededor de la frecuencia de Nyquist se elimina por completo, y la figura 7 muestra la respuesta en frecuencia de la respuesta impulsional de la figura 5.

• 
  Figura 6. Respuesta en frecuencia con frecuencia digital antes del filtrado.

• 
  Figura 7. Respuesta en frecuencia con frecuencia digital después del filtrado.

• Wang, Yanghua (2008). Filtrado Q inverso sísmico . Blackwell Pub. ISBN 978-1-4051-8540-0.

• Filtrado Q inverso estabilizado por Knut Sørsdal