Facultat de Ciències
Memòria del Treball de Fi de Grau
Estudio de oscilaciones transversales de espículas cromosféricas empleando “Empirical
Mode Decomposition”
Matheus Aguiar-Kriginsky Silva Grau de Física
Any acadèmic 2017-18
Treball tutelat per Ramón Oliver Herrero Departament de Física
S'autoritza la Universitat a incloure aquest treball en el Repositori Institucional per a la seva consulta en accés obert i difusió en línia, amb finalitats exclusivament acadèmiques i d'investigació
Autor Tutor
Sí No Sí No
X X
Paraules clau del treball:
Empirical Mode Decomposition, Física solar, Métodos estadísticos para señales multicomponente, Oscilaciones transversales de espículas.
Índice
1. Introducción 4
1.1. La Atmósfera solar . . . 4
1.1.1. La fotosfera . . . 4
1.1.2. La Cromosfera. . . 4
1.1.3. La Corona solar . . . 5
1.2. Las Espículas Solares . . . 5
1.3. Objetivo del trabajo. . . 6
2. Fast Fourier Transform (FFT) 7 3. Descomposición Empírica de Modos 8 3.1. Funciones Intrínsecas de Modo. . . 8
3.2. Algoritmo EMD . . . 9
3.3. Ensemble Empirical Mode Decomposition (EEMD). . . 13
3.4. Complete Ensemble Empirical Mode Decomposition With Adaptive Noise (CEEMDAN) . . . 14
4. Resultados 15 4.1. Análisis con la FFT . . . 16
4.2. Análisis con EMD . . . 18
4.3. Análisis con EEMD . . . 21
4.4. Análisis CEEMDAN . . . 24
4.5. Residuo y Tendencia . . . 26
5. Conclusiones y Observaciones 28
Referencias 29
1. Introducción
A parte de establecer posiblemente la primera clasificación estelar y producir un mapa exac- to del cráter lunar Copernicus, en el año 1877 Angelo Secchi (29 de Junio de 1818 - 26 de Febrero de 1878) publicóLe Stelle(Las Estrellas), donde se demuestra por primera vez que las protuberancias observadas durante los eclipses solares pertenecen efectivamente al Sol, y se relata la presencia de manchas sobre su superficie. Alrededor de estas manchas, Secchi relata como se pueden observar pequeños chorros que emergen desde capas más profundas del Sol, de corta duración y de dinámica variable.
Desde entonces, estos chorros (denominados espículas) se siguen observando y estudiando como parte de la dinámica solar.
Las espículas son chorros de plasma que se elevan a velocidades de 20 km/s, pudiendo al- canzar alturas de miles de kilómetros hasta desaparecer. Como ya apuntaba Secchi, son de corta duración, típicamente entre 5 y 15 minutos.
1.1. La Atmósfera solar
La atmósfera solar se suele subdividir en tres capas, de diferente grosor y con características específicas. Para tener una visión global de cómo está dividida, se presenta un breve resu- men.
1.1.1. La fotosfera
La fotosfera es la capa inferior, con unos 300 - 500 km de grosor. Es la capa más fría, con un descenso de la temperatura con la altura. En esta capa es donde se pueden observar las fáculas solares, los gránulos y las manchas solares (ver la figura1.
Figura 1: Imagen de la fotosfera solar (1992), la capa más conocida. Se puede apreciar la presencia de man- chas solares. Verhttps://solarscience.msfc.nasa.gov/surface.shtml.
1.1.2. La Cromosfera
Es la capa intermedia, con un grosor de entre 3000 y 5000 km. Al contrario de lo que ocurre en la fotosfera, la temperatura de la cromosfera aumenta con la altura, desde unos 6000 K en la base hasta unos 20000 K en la parte superior. Estas temperaturas propician la emisión Hαdel hidrógeno, que se caracteriza por un color rojizo, y esto es lo que le da su nombre a la capa (cromosfera = esfera de color). Por encima de la cromosfera se pueden observar estructuras cuyas propiedades (densidad y temperatura) son típicamente cromosféricas: las espículas (que serán descritas en más detalle en la sección1.2) y las protuberancias, que se proyectan sobre el limbo solar durante eclipses totales del Sol (figura2) o sobre el disco solar.
1.1.3. La Corona solar
Consiste en la capa exterior de la atmósfera solar, extendiéndose millones de kilómetros en el espacio (figura2). Es visible durante los eclipses totales o con un coronógrafo, que con- siste en un telescopio dotado de un disco para tapar la fotosfera y producir un eclipse Solar artificial. Las estructuras más destacadas que se observan en la corona solar son los denomi- nados agujeros coronales, que son zonas de menor densidad y temperatura en comparación al resto de la corona, así como los bucles coronales.
Figura 2: Imagen del Sol durante el eclipse solar del 11 de Julio de 1991. La emisiónHαque cau- sa el color rojizo se puede apreciar sobre el limbo solar y proviene de las protuberancias. Las es- tructuras blancas que se extienden hacia el espacio exterior constituyen la corona solar. Ver https://solarscience.msfc.nasa.gov/chromos.shtml.
1.2. Las Espículas Solares
Las espículas son observadas constantemente en la cromosfera solar como material que es lanzado desde la parte baja de la atmósfera solar hacia arriba, con una presencia de unas 300.000 espículas en todo el Sol en cualquier momento. La figura3muestra que las espículas proyectadas sobre el limbo solar destacan en imágenes tomadas en diferentes longitudes de onda alrededor del centro de la línea Hα, mientras que las que se proyectan sobre el disco son visibles únicamente para ciertas longitudes de onda en las alas de esta línea espectral.
Análisis detallados de las espículas revelan que existen al menos dos tipos, con caracterís- ticas diferentes. Las espículas de tipo I tienen dinámicas dominadas por escalas de tiempo de entre 3 y 7 minutos, presentando un movimiento de subida (observable en esteenlace [6]) y de bajada (apreciable en esteenlace[6]). Son las espículas que más abundan en las regiones activas del Sol, abarcando hasta el 60 % de las espículas observadas en algunos es- tudios. Como tubos de plasma, tienen un diámetro de entre 400 y 1500 km, y una longitud de 5000-10000 km. Las denominadas espículas de tipo II se caracterizan por no tener una fase de descenso, con unas escalas de tiempo de menos de 3 minutos. También suelen tener un diámetro menor (≤200 km), y a bajas alturas parecen ser las que dominan.
El origen de las espículas solares ha sido un tema de debate desde su descubrimiento. En un estudio reciente, Martínez-Sykora et al. (2017) [5] proponen que la clave está en la pre- sencia de partículas neutras en el plasma, lo que es posible debido a que la temperatura es insuficiente para tener el gas completamente ionizado.
El modelo propuesto se basa en la presencia de convección en la zona por debajo de la su- perficie solar, donde se generan puntos neutros del campo magnético (donde el campo se anula). La convección eleva dichos puntos hasta la superficie, y una vez allí el punto se des- hace, liberando energía y acelerando el plasma en el proceso. Las espículas son chorros de
Figura 3: Serie de imágenes hechas en 15 longitudes de onda alrededor del centro de la líneaHα; las 7 pri- meras corresponden al ala azul de la línea, la octava al centro de línea y las 7 siguientes al ala roja. Se pue- den apreciar las espículas como las estructuras alargadas, observadas en emisión sobre el limbo solar y en absorción sobre el disco solar. Se puede ver la película formada por estas imágenes en esteenlace. Ver https://arxiv.org/pdf/1602.05185.pdf.
Sin las partículas neutras en el modelo (incluirlas eleva el coste computacional considera- blemente), no susceptibles al campo magnético, los puntos no se elevan hasta la superficie solar. Además, la fricción entre partículas neutras y los iones calienta más el plasma dentro y alrededor de las espículas, lo que permite elevar su temperatura hasta los valores observa- dos.
Martínez-Sykora et al. (2017) [5] compararon los resultados de sus simulaciones con obser- vaciones del Interface Region Imaging Spectrograph (IRIS), resultando en concordancias sa- tisfactorias.
Aunque su papel en la dinámica solar puede ser distinto, tanto las espículas de tipo I como II presentan una característica común: oscilan transversalmente. Dicha oscilación lleva siendo objeto de estudio desde su descubrimiento. Y parece ser fruto de la propagación de ondas a lo largo de las espículas, en modos de oscilación conocidos como "kink waves". Dada su geometría cilíndrica y superior densidad relativa al medio que la rodea, una espícula puede servir de guía de ondas para la propagación de perturbaciones generadas en la fotosfera.
La fuente de estas perturbaciones, al menos en el caso de las espículas de tipo I, podría ser la fuga de ondas de presión (o denominadas modo p) producidas en el interior solar. El es- tudio del período de oscilación transversal de las espículas por lo tanto deviene crucial para entender su naturaleza.
1.3. Objetivo del trabajo
El objetivo de este trabajo es estudiar las oscilaciones transversales de las espículas solares introduciendo un método alternativo al análisis de Fourier. El método, denomidado Des- composición en Modos Empíricos (EMD, del inglésEmpirical Mode Decomposition), consis- te en descomponer una función en un conjunto de Funciones Intrínsecas de Modo (IMFs, del InglésIntrinsic Mode Functions). Estos modos forman no solo una base de autovecto- res mucho menor que la correspondiente descomposición en funciones trigonométricas del análisis de Fourier, sino que cada uno de los autovectores tiene un significado más simple de interpretarse, lo que hace el análisis más sencillo.
2. Fast Fourier Transform (FFT)
La descomposición de una función periódica en una suma de funciones trigonométricas se lleva haciendo desde la invención de la transformada de Fourier discreta para saber qué frecuencias están presentes en una señal. Esta descomposición se basa en el hecho de que las funciones trigonométricas (o, en general, las exponenciales complejas) forman una base ortogonal y completa en el espacio de funciones en el que viven, y por lo tanto pueden repro- ducir idénticamente la función original si se dispone de todos los coeficientes del desarrollo.
La transformada de Fourier de una función f(t) se define como:
fˆ(k)= Z ∞
−∞
f(t)e−2πi ktd t. (1)
Al hacer la transformada de Fourier, se obtiene una función en el espacio de frecuencias. Su operación inversa es la antitransformada, que lleva a la función original f(t):
f(t)= Z ∞
−∞
fˆ(k)e2πi ktd k. (2)
Si la función f(t) es periódica, la transformada de Fourier no es más que una serie de deltas de Dirac con picos en las frecuencias propias de la señal. Si la señal no es continua, sino que se dispone de un número finito de puntos equiespaciados temporalmente (tk=k∆, con k=0, ...,N−1) se realiza la transformada de Fourier discreta (DFT):
Fn=
N−1
X
k=0
fke−2πi nk/N. (3)
En general la DFT de una secuencia deN números reales es una secuencia deN números complejos.
La denominada Fast Fourier Transform (FFT) es un algoritmo para calcular la DFT reducien- do el número de computaciones necesarias para transformarNpuntos de 2N2a 2NlnN. De una manera u otra, algoritmos similares al de la FFT se utilizan desde 1805 con Gauss (H.
H. Goldstine [2] le atribuye un algoritmo similar al de la FFT para calcular los coeficientes de una serie de Fourier finita). Sin embargo, el análisis espectral de Fourier tiene grandes des- ventajas: el sistema descrito ha de ser lineal, y los datos han de ser periódicos o estacionarios.
Si no se dan estas circunstancias, el espectro de frecuencias resultante tiene poco sentido fí- sico. A parte de esto, al hacer la descomposición de la señal en funciones periódicas:
F(t)=Re à N
X
n=−N
cNe2πi nt
!
, (4)
puede que sean necesarios demasiados términos en la suma para poder reproducir de ma- nera satisfactoria la señal original, algo que puede ser costoso y poco práctico. Esta depen- dencia con el número de términos que se elige en la suma se puede apreciar en la figura4. Se muestra una función compuesta de cuatro escalones, formada por 350 puntos. Claramente se necesitan una gran cantidad términos para obtener una descomposición que se ajuste bien a la señalF(t).
Figura 4: Suma de la serie de Fourier truncada a diferentes valores deNde una función formada por cuatro escalones (línea naranja) comparada con la función original (línea azul).
3. Descomposición Empírica de Modos
3.1. Funciones Intrínsecas de Modo
Dadas las dificultades presentadas por la descomposición habitual en series de Fourier de una señal con componentes oscilatorias no estacionarias, en este proyecto se propone el uso de un nuevo método para el tratamiento de señales (Huang et al. 1998 [3]). Se basa en descomponer la serie de datos en las denominadas Funciones Intrínsecas de Modo (IMF, del inglésIntrinsec Mode Functions).
Una IMF es cualquier función que cumple las siguientes condiciones:
1. El número de extremos y el número de veces que cruza el 0 ha de ser el mismo o diferir en uno.
2. En cualquier punto, el valor medio de la envolvente definida por los máximos locales y la envolvente definida por los mínimos locales es 0.
Como se puede observar, las funciones base de la DFT (o sea, los senos y cosenos con una frecuencia fija) son un caso particular de IMF, aunque la definición dada incluye funciones mucho más generales, cuya frecuencia y amplitud pueden variar con el tiempo.
La motivación del método de Descomposición Empírica de Modos (EMD, del InglésEmpi- rical Mode Decomposition) es poder tratar con sistemas no lineales y no estacionarios y, por lo tanto, con señales cuya frecuencia instantánea depende del tiempo.
Dada una IMFX(t), su frecuencia instantánea se calcula de la siguiente manera:
1. Se realiza la transformada de Hilbert de la IMF:
Y(t)= 1 πV.P.
Z +∞
−∞
X(t0)
t−t0d t0 (5)
donde V.P. indica el valor principal de Cauchy. Esta definición resulta en una función Y(t) tal que la función
Z(t)=X(t)+i Y(t)=a(t)eiθ(t) donde
a(t)=p
X2(t)+Y2(t), θ(t)=arctan
³Y(t) X(t)
´ , es una función analítica.
2. Se define la frecuencia instantánea como:
ν(t)= 1 2π
dθ(t)
d t (6)
La condición 2 ideal para la definición de las IMFs sería que ‘la media local de la función sea cero’, pero dicha definición deja de tener sentido si la señal no es estacionaria. Esta aproxi- mación en la definición de las IMFs puede introducir errores en la frecuencia instantánea para señales no lineales, pero aun así la frecuencia instantánea será consistente con la física del sistema en estudio.
Definidas las IMFs y la frecuencia instantánea, una señal puede ser descompuesta en una base casi ortogonal (la ortogonalidad es una aproximación pero bastante cercana a la reali- dad) de IMFs, cada una con una frecuencia instantánea característica que viene influencia- da por la física del problema y por el posible ruido. Se obtiene además un residuo que ya no puede ser descompuesto en IMFs.
Ahora se trata de definir un algoritmo para encontrar las IMFs propias de cada señal que se vaya a analizar. A continuación se presentan 3 algoritmos, cada uno con mejoras respecto a los anteriores.
3.2. Algoritmo EMD
Este proceso para la obtención de las IMFs propias en la señal parte de las siguientes premi- sas [3]:
La señal tiene amenos dos extremos, un máximo y un mínimo.
La escala de tiempo característica viene definida por el intervalo de tiempo entre dos extremos sucesivos.
Si la señal no tuviera máximos ni mínimos pero sí puntos de inflexión, se pueden ob- tener los extremos derivando y luego integrando los resultados obtenidos.
Las condiciones anteriores resumen la metodología a seguir, ya que si hay muchos modos de oscilación superpuestos la manera de extraerlos es fijándose en la escala de tiempo que hay entre dos extremos consecutivos del modo de oscilación respectivo. Es una forma bastante intuitiva de proceder, ya que al observar una serie temporal lo primero que uno observa para saber si hay diferentes modos de oscilación es cómo van variando los extremos.
Para explicar el funcionamiento del algoritmo de la EMD, se utilizará la función siguiente:
f(t)=cos(2π0,5t)+cos(2π0,05t). (7) La función7se representa en la figura5, y está claro por inspección de los máximos y míni- mos que hay dos modos de oscilación presentes.
Figura 5: Señal temporal de la función7.
Teniendo presente la definición de las IMFs, el primer paso es identificar la envolvente de- finida por los máximos de los datos y la envolvente definida por los mínimos. Se ajusta un
‘spline’ cúbico a las dos envolventes, y punto a punto se calcula la media entre los dos ajustes (m1). Se resta la media de los datos:
f(t)−m1(t)=A1(t) (8)
En este caso, con una iteración ya se llega a una IMF (Figura6). Sin embargo, debido a po- sibles problemas con el ajuste de los ‘splines’ (pueden desarrollar colas en los extremos) o a la posible generación de nuevos extremos al restar la media puede ser necesario repetir el proceso iterativamente hasta que se obtenga una IMF de verdad. A A1se la trataría de la misma manera que f(t), ajustando ‘splines’ cúbicos y restando la media punto a punto
A1(t)−m11(t)=A11(t) (9)
y se comprobaría si A11(t) es o no una IMF.
Figura 6: Ajuste de ‘splines’ cúbicos sobre los extremos de los datos (izquierda) y el resultado de restar la me- dia a los datos originales (derecha). Se obtiene una IMF en la primera iteración.
Obtenida la primera IMF después deniteraciones,
A1(n−1)(t)−m1n(t)=A1n, (10)
se define
c1(t)=A1n(t) (11)
como la primera IMF de los datos.
Aunque este método funciona para extraer las IMFs, si se hacen demasiadas iteraciones se pueden eliminar posibles fluctuaciones de amplitud y de frecuencia que son característi- cas de la física de la situación. Por lo tanto, el proceso se ha de parar con algún criterio. El propuesto por Huang et al. (1998) [3] es limitar el tamaño de la desviación estándar (DE) calculada entre dos pasos consecutivos:
DE=
T
X
t=0
"¯
¯A1(n−1)(t)−A1n(t)¯
¯
2
A21(n−1)(t)
#
(12)
Los valores típicos de DE propuestos son de 0.2 - 0.3, aunque conviene mirar caso a caso.
Siguiendo con el proceso, una vez obtenida la IMF c1 (de menor período) se resta de los datos:
f(t)−c1(t)=r1(t) (13)
Conr1se procede igual que anteriormente para obtenerc2, y así sucesivamente hasta obte- ner un último residuo:
f(t)−
n
X
k=1
ck(t)=rn(t) (14)
El residuornya es una función monotónica de la cual no se pueden extraer más IMFs, o es de valor demasiado pequeño en comparación a un valor predeterminado. Se obtiene entonces una descomposición de la señal enn modos empíricos (es decir,n IMFs) y un residuo que es una posible tendencia de los datos o una constante:
f(t)=
n
X
k=1
ck(t)+rn(t). (15)
El número de IMFs en que se descompone una señal (n) es mayor cuanto más elevado es el número de muestras de la señal y depende del ruido contenido en la señal, del número de modos físicos que intervienen, del valor de DE adoptado, etc. En el caso del ejemplo pro- puesto en la figura5, se obtienen dos IMFs y un residuo. La primera IMF se muestra en el panel de la derecha de la figura6, mientras que la segunda IMF y el residuo se presentan en la figura7.
Figura 7: Segunda IMF y el residuo del análisis propuesto como ejemplo.
Para completar el análisis y entender el sentido físico de cada IMF, en la figura8se presentan los resultados del cálculo de la frecuencia instantánea de la ecuación6. En el panel de la izquierda se muestra ν(t) para todo el rango temporal, y se puede apreciar los efectos de frontera causados por la transformada de Hilbert. Estos efectos no son físicos, y por eso en el panel de la derecha se hace una ampliación de los resultados a la zona central. Al analizar la zona central, se observa como la frecuencia instantánea de cada IMF es la correspondiente a la frecuencia de oscilación de la función de la que originan los datos (ecuación7), y por lo tanto es intrínseca a la naturaleza de la señal.
Figura 8: Izquierda:ν(t) para las dos IMFs del ejemplo propuesto. Derecha: ampliación de la imagen de la izquierda sobre la zona alejada de las fronteras
El algoritmo descrito funciona bien para el ejemplo propuesto, pero puede generar proble- mas en casos más complejos. El mayor inconveniente es la posibilidad de la mezcla de mo- dos. Este fenómeno se da cuando dos oscilaciones con diferentes escalas de tiempo están presentes en una misma IMF, o cuando oscilaciones con la misma escala de tiempo están presentes en dos IMFs diferentes. Esto se debe a la naturaleza local de la EMD, ya que el primer modo contiene la frecuencia localmente más alta y así sucesivamente hasta llegar a los modos de frecuencia más baja. Un ejemplo claro de mezcla de modos se puede apreciar en la figura9, en la que se presenta una señal compuesta por un tono puro sostenido y otro intermitente. La primera IMF recoge la frecuencia localmente más alta, que cambia con el tiempo debido a la señal descompuesta. Esto hace que en dos IMFs diferentes esté presente la misma frecuencia instantánea, como se puede apreciar en el panel inferior central de la figura9.
Figura 9: Representación de la señal correspondiente a un tono puro sostenido y otro intermitente (arriba, panel izquierdo), así como las dos IMFs resultantes de la descomposición EMD (arriba, paneles central y de- recho). También se representan los diagramas frecuencia-tiempo para la señal (abajo, panel izquierdo) y para las IMFs (abajo, paneles central y derecho) que evidencian la mezcla de modos en la primera IMF (etiquetado modo 1). Fuente:[1].
Otra característica particular del método EMD es que los autovectores que constituyen la ba- se en la que se descompone la señal no están prefijados, sino que vienen determinados por la propia señal. Esto lo diferencia de otros métodos como la DFT o la CWT (transformación continua de onditas, del inglésContinuous Wavelet Transform).
3.3. Ensemble Empirical Mode Decomposition (EEMD)
Para poder solucionar el problema de mezcla de modos presente en el algoritmo original, una propuesta (Wu, et al. (2009)[8]) es utilizar ruido generado aleatoriamente para ayudar a separar mejor las escalas de tiempo que en el método EMD. Se añade ruido blanco a la señal diversas veces, y como el ruido es cada vez diferente las IMFs de una realización no son las mismas que las IMFs de realizaciones con diferentes ruidos. Con un número adecuado de intentos realizados, el ruido añadido se puede eliminar haciendo una media del colectivo de IMFs obtenidas en las diferentes realizaciones. El procedimiento es el siguiente:
1. A partir de la señal original f(t) se generanNrealizaciones de ruido blanco aleatorio:
Fi=f(t)+wi (i=0,...,N), dondewi es cada realización de ruido.
2. Se descompone cadaFi en IMFs con el método EMD descrito anteriormente.
3. Las diferentes IMFs obtenidas se promedian sobre todo el conjunto de descomposi- ciones realizadas, y los promedios se toman como las IMFs verdaderas:
I M Fk= 1 N
N
X
i=1
I M Fki (16)
dondekrepresenta la IMF a ser promediada de lasK obtenidas (k=1,2,...,K)
La descomposición se calibra fijando la amplitud del ruido añadido y el tamaño del colectivo de realizaciones (N) para obtener los mejores resultados posibles.
Aunque este método puede ayudar a solucionar el problema de la mezcla de modos, parece evidente que añadir ruido a la señal y hacer la media puede dejar ruido residual sobre las IMFs. Otro problema es que diferentes realizaciones de la descomposición EMD sobre dife- rentesFipueden dar lugar a un número diferente de IMFs, por lo cual el promedio realizado se vería comprometido.
3.4. Complete Ensemble Empirical Mode Decomposition With Adaptive Noise (CEEMDAN)
Parece necesaria una mejora del método EEMD presentado para que no se den los errores que introduce al intentar retirar el problema de la mezcla de modos. Con ello se presenta una modificación del método (Torres, et al.(2011) [7]). En el algoritmo EEMD, cada realizaciónFi se descompone independientemente de las demás, resultando en residuos diferentes para cada realización:
rki =rki−1−I M Fki (17)
La propuesta del método CEEMDAN es obtener un primer residuo:
r1=f(t)−I M F1 (18)
dondeI M F1se calcula con el algoritmo EEMD presentado antes. Obtenidor1, se calcula el primer modo EMD (la primera IMF de la descomposición) sobre un colectivo de realizacio- nes der1más ruido, de manera algo similar a la del método EEMD, y al promediar se obtiene I M F2. DefiniendoEj(·) como un operador que actuando sobre una señal produce el modo j obtenido con EMD, el algoritmo general es el siguiente:
1. Se descomponenNrealizaciones de f(t)+ε0wi (i=1,...,N) con el algoritmo EMD para calcular:
I M F1= 1 N
N
X
n=1
I M F1n (19)
2. Se calcula el primer residuo:
r1=f(t)−I M F1 (20)
3. Se descomponenN realizaciones der1+ε1wi y la segunda IMF se define a partir del primer modo obtenido haciendo la descomposición EMD:
I M F2= 1 N
N
X
n=1
E1(r1+ε1wi) (21)
4. Se obtiene el residuor2=r1−I M F2
5. Se procede de la misma manera hasta que del futuro residuo no se puedan extraer más IMFs.
Los coeficientesεi se utilizan para controlar la relación señal a ruido, y los valores de esta amplitud se fijan de acuerdo con la frecuencia de la señal a ser descompuesta.
4. Resultados
Para el estudio de la oscilación de las espículas en la cromosfera solar, en este proyecto se utilizarán los datos obtenidos por el instrumento Rapid Oscillations in the Solar Atmosp- here (ROSA) presente en el Dunn Solar Telescope (National Solar Observatoy, NM, EEUU).
Consisten en una serie de observaciones hechas entre las 13:46 y 14:40 UT el 28 de Mayo de 2009 (una animación se puede ver en esteenlace), con una cadencia efectiva de 4.2243 s. Los datos originan de Kuridze et al. (2013) [4], que trataron los datos y luego hicieron un análisis con el cuál se pueden hacer comparaciones.
Para estudiar la oscilación transversal de una espícula, se ha de elegir una que esté más o menos aislada dentro de todas las que se observan (figura10). Una vez elegida la espícula a estudiar se realizan una serie de cortes horizontales sobre esta (figura10).
Figura 10: Izquierda: Vista de la imagen hecha a partir de los datos obtenidos con ROSA, donde se recuadra la zona en la que se identifica una espícula útil para el análisis. Derecha: Ampliación del recuadro rojo del panel izquierdo, así como de 5 de los cortes hechos para analizar la oscilación de la espícula.
Para un instante de tiempo fijo, un corte cualquiera permite apreciar la variación de la in- tensidad en la dirección transversal a la espícula en un punto de esta. Si se apilan las series temporales de la intensidad frente a la posición para un corte concreto, entonces se obtienen los resultados de la figura11, izquierda. Se observa claramente el desplazamiento transver- sal de la espícula en cada uno de los cinco puntos. A continuación se considera que, para cada instante de tiempo, el mínimo de intensidad en un corte corresponde a la posición de la espícula. La posición del mínimo se determina haciendo un ajuste gausiano, que permi- te mejor resolución espacial de la que se dispone solo con el corte. La representación del mínimo de intensidad frente al tiempo se presenta en la figura11, derecha, para los cinco cortes.
Figura 11: Izquierda: Evolución de la intensidad de los cortes señalados en el panel derecho de la figura10.
Los ejes horizontal y vertical corresponden al tiempo (t) y la posición transversal respecto de la espícula (s).
Derecha: resultado del ajuste realizado a cada paso de tiempo sobre los cortes.
La espícula elegida se ha observado durante 110 pasos de tiempo, que equivalen a unos 7.8 minutos de oscilación (observable en esteenlace, en el que la oscilación transversal se apre- cia a simple vista). Las series temporales resultantes presentan una oscilación completa, lo que hace de la espícula una buena candidata a ser tratada con el análisis de Fourier y tam- bién con los métodos presentados anteriormente (EMD, EEMD, CEEMDAN).
4.1. Análisis con la FFT
Aunque en la figura10se muestran 5 series, en la realidad se han hecho 10 cortes (no se han pintado los demás por claridad). El tratamiento con la transformada de Fourier para saber cuál el período de cada una de las series consiste en aplicar la FFT a cada una de ellas y ver qué frecuencia es la que más contribuye.
La resolución en frecuencia es de 0.00215 Hz, lo que equivaldría a un período de unos 460 segundos. Por lo tanto, como el período que se intenta medir es del orden de 200 segundos (como se aprecia en las series temporales de los diferentes cortes), hay problemas de resolu- ción. El error en la determinación de la frecuencia de una oscilación se puede estimar como el espaciado entre frecuencias, es decir:∆ν=0,00215 Hz. Este error sobre la frecuencia se propaga al valor de la estimación del período, que tendría un error∆P=P2∆ν. Para perío- dos del orden de 200 segundos,∆P/Pes de un 43 %. Esto quiere decir que la indeterminación sobre el valor del período es elevada. La transformada de Fourier de uno de los cortes se pue- de ver en la figura12. Haciendo la media de todas las frecuencias, se obtiene un período de 201segundos.
Figura 12: Representación de la transformada de Fourier de la serie temporal de uno de los cortes.
Para observar la necesidad de muchos modos para describir la señal en términos de una serie de Fourier, se muestra en la figura13la representación en series de Fourier de uno de los cortes. La necesidad de muchos términos puede ser costosa y además cada modo da poca información acerca de la física del sistema siendo estudiado.
Figura 13: Suma de la serie de Fourier truncada a diferentes valores deNde la serie temporal de uno de los cortes (línea amarilla) comparada con la función original (línea azul). Solo cuandoN =110 (el número de puntos de la serie temporal), se obtiene una recostrucción correcta de la serie.
4.2. Análisis con EMD
Para comprobar la utilidad del Empirical Mode Decomposition, se analizan las 10 series tem- porales, fijandoDE=0,1. A modo de ejemplo, una de las descomposiciones, la del primer corte de la figura10, se puede ver en la figura14.
Figura 14: Izquierda: Serie temporal de uno de los cortes mostrados en la figura11. Derecha: Descomposición de la serie temporal en IMFs con el algoritmo EMD.
Con el algoritmo EMD, la señal temporal ha sido descompuesta en 5 IMFs (la sexta compo- nente en la figura14es el residuo). Para comprobar que la descomposición es completa, en la figura15se compara la señal a la suma de las IMFs y el residuo, y se puede ver como la suma y la serie son idénticas, con una máxima discrepancia del orden de 10−15.
Figura 15: Suma de las IMFs (con el método EMD) y señal temporal comparadas.
Con el nuevo método descrito, se ha descompuesto una señal en muchas menos autofun- ciones que con la serie de Fourier, y la señal se reproduce idénticamente. Otra de las grandes ventajas que presenta este método es que las IMFs en las que se descompone la señal tem- poral tienen un sentido físico mucho más fácil de identificar que en el caso del análisis de Fourier. El objetivo del algoritmo EMD en este proyecto es caracterizar el período de oscila- ción de las espículas, y si se compara la señal temporal en la figura14con sus IMFs, se puede apreciar que la IMF número 4 tiene un aspecto muy similar al de la señal, es decir, representa el modo de oscilación intrínseco de la espícula.
Ambas funciones se comparan en la figura16. Salvo la tendencia que presenta la espícula (no presente en la IMF pero absorbida en el residuo), se aprecia la similitud entre ambas funciones.
Figura 16: Izquierda: Serie temporal de uno de los cortes mostrados en la figura10. Derecha: IMF 4 obtenida con el método EMD. La semejanza es evidente, salvo la tendencia temporal y la variabilidad de alta frecuen- cia.
Debido a problemas que presenta la transformada de Hilbert en las fronteras, se observan fluctuaciones deν(t) que no son reales hacia las fronteras de los datos. Si el objetivo es obte- ner el período de la oscilación de la espícula, lo ideal es fijarse en los valores más centrales de la frecuencia, debido a los problemas de fronteras mencionados. En el panel izquierdo de la figura17se presentan las frecuencias instantáneas de las IMFs que representan el modo fí- sico de oscilación de la espícula para las 10 series temporales estudiadas, mientras que en el panel derecho de la misma figura se hace una ampliación a la zona más central. Los efectos de frontera son visibles, provocando oscilaciones que no son propias del sistema observado.
Sin embargo, en la zona central, en todas las series temporales el valor de la frecuencia ins- tantánea parece ser el mismo, de aproximadamente 0.005s−1, lo que equivale a un período de 200 segundos.
Figura 17: Izquierda:ν(t) para la IMF (obtenida con el método EMD) que caracteriza la oscilación transversal de la espícula para todas las series temporales estudiadas. Derecha: Ampliación de las figuras de la izquierda en el intervalo de tiempo 100≤t≤250s.
Para comprobar que utilizando solamente la descomposición con el algoritmo EMD se pue- de producir mezcla de modos, en la figura18 se representan las frecuencias instantáneas de las IMFs 3, 4 y 5 de la serie temporal descrita en la figura 16. Se hace una ampliación de la zona central ya que la mezcla en los extremos puede deberse a los problemas de la transformada de Hilbert y no al método utilizado. Se puede observar como hay una misma frecuencia presente en dos IMFs diferentes, siendo una de ellas la que se ha utilizado para caracterizar la oscilación de la espícula. Esto podría afectar a los resultados obtenidos, y sirve de motivación para la utilización de los otros métodos descritos.
Figura 18: Frecuencia instantánea de las IMFs 3, 4 y 5 obtenidas con el método EMD. Hay un instante de tiem- po en el que dos IMFs tienen la misma frecuencia.
4.3. Análisis con EEMD
Se tratan los datos obtenidos de la misma manera que al utilizar el algoritmo EMD, solo que a la hora de hacer la descomposición en IMFs se utiliza el algoritmo EEMD. Se utiliza un colectivo deN=500 realizaciones, y con un valorε= 0.02. La descomposición de la misma señal temporal que la de la figura 14se muestra en la figura 19. Como se puede observar, se obtienen más IMFs que con el método anterior. Esto se debe al ruido añadido para ha- cer la descomposición, y no tiene por qué significar que hay o no más modos de oscilación físicamente significativos que el algoritmo EMD no es capaz de extraer. La completitud del método EEMD se puede comprobar en la figura20.
Figura 19: Izquierda: Serie temporal de uno del panel izquierdo de la figura16. Derecha: Descomposición de la serie temporal en IMFs con el algoritmo EEMD.
Figura 20: Suma de las IMFs (con el método EEMD) y señal temporal comparadas.
Análogamente al análisis hecho en el apartado anterior, se identifica la IMF 4 como la que refleja la oscilación de la espícula (figura21). También se compara la IMF 4 obtenida con EEMD con la IMF 4 obtenida en la misma serie temporal con el método EMD. Al comparar- las, parece que la adición de ruido ha afectado la amplitud de la IMF, como se puede ver al principio.
En la figura22se muestraν(t) para las IMFs de los 10 cortes que corresponden con el modo de oscilación de la espícula (panel izquierdo), así como una ampliación a la zona central en el panel izquierdo. Los valores deν(t) también parecen centrarse en el valor 0.005s−1, pero la oscilación respecto de este valor es superior al caso del análisis con EMD (figura17).
Figura 21: Izquierda: Serie temporal del panel izquierdo de la figura19. Derecha: IMF 4 obtenida con EEMD (en azul), comparada con la IMF 4 obtenida con EMD (en naranja).
Figura 22: Izquierda:ν(t) para la IMF (obtenida con el método EEMD) que caracteriza la oscilación de la es- pícula para todas las series temporales estudiadas. Derecha: Ampliación de las figuras de la izquierda sobre la zona central.
Para ver si el método EEMD resuelve la mezcla de modos del corte estudiado presente en la descomposición hecha con el método EMD, se representan las ν(t) de las IMFs 3, 4 y 5 de la figura19en la zona central, como se hizo en el caso del análisis con el método EMD (figura 18). Se puede apreciar (figura23) que la mezcla de modos en la zona donde había anteriormente ya no se produce, pero debido a la introducción del ruido blanco se mezclan los modos en otro instante.
Figura 23: Comparación deν(t) para las IMFs 3, 4 y 5 obtenidas con EEMD.
4.4. Análisis CEEMDAN
Se ha constatado que el algoritmo EMD presenta mezcla de modos, y que aunque el método EEMD pueda solucionar dicha mezcla, pueden aparecer nuevas mezclas y además se obtie- ne un número distinto de IMFs. Para comparar, se vuelve a descomponer la señal temporal del mismo corte en IMFs con el método CEEMDAN (figura24). Se ha utilizado un colectivo de 500 realizaciones, con un valor deεconstante de 0.002. En comparación con los métodos anteriores, se obtienen 3 IMFs y un residuo. Esto puede significar (o no) que los métodos anteriores extraigan modos de oscilación que no tienen sentido físico.
Figura 24: Izquierda: Serie temporal de uno del panel izquierdo de la figura16. Derecha: Descomposición de la serie temporal en IMFs con el algoritmo CEEMDAN.
Comparando la señal temporal con las IMFs de la descomposición, la que parece correspon- der al modo de oscilación de la espícula es la IMF 3 (figura25). Al comparar esta IMF con las correspondientes a los otros métodos (figura25, derecha), se puede observar que utilizar el método CEEMDAN no afecta casi a la amplitud de la IMF, y esta es más similar a la original (obtenida con el método EMD).
En la figura26se representa la frecuencia instantánea de las IMFs obtenidas con el método CEEMDAN correspondientes al modo de oscilación de la espícula para las 10 series tempo- rales, y oscilan alrededor del mismo valor que con los otros métodos, pero sin oscilar tanto como con el método EEMD.
Figura 25: Izquierda: Serie temporal. Derecha: IMF 3 obtenida con CEEMDAN (en azul), comparada con la IMF 4 de la figura16obtenida con EMD (naranja), y la IMF 4 de la figura21obtenida con EEMD (verde).
Figura 26: Izquierda:ν(t) para la IMF (obtenida con el método CEEMDAN) que caracteriza la oscilación de la espícula para todas las series temporales estudiadas. Derecha: Ampliación de las figuras de la izquierda sobre la zona central.
EMD y EEMD en la región central, se representaν(t) para las IMFs 2 y 3 (figura27). Como se puede apreciar, la mezcla de modos no solo no está presente, sino que las frecuencias se alejan en todo el dominio de interés, lo que significa que el modo de oscilación que una representa no coinciden con el de la otra, una de las premisas de la descomposición en IMFs.
Figura 27: Comparación deν(t) para las IMFs 3 y 4 de la figura25obtenidas con CEEMDAN.
La completitud del método CEEMDAN se puede comprobar en la figura28.
Figura 28: Suma de las IMFs (con el método CEEMDAN) y señal temporal comparadas.
.
4.5. Residuo y Tendencia
En el análisis realizado en el artículo del que provienen los datos (Kuridze et al. [4]), estos au- tores restan una tendencia lineal a cada serie temporal antes de realizar el análisis de Fourier sobre los datos. Con ello, obtienen un período de oscilación medio de 180 segundos, que no parece estar alejado de lo que se percibe a simple vista. Sin embargo, esta tendencia lineal que restan puede afectar a sus resultados. En este trabajo se ha procedido de otra manera, utilizando los métodos EMD, EEMD y CEEMDAN. Para utilizar dichos métodos, no es nece- sario manipular los datos restando ninguna tendencia, ya que la tendencia que puede estar presente en la señal temporal resulta en el residuo de cada descomposición. Además, las ten- dencias resultantes de aplicar los tres métodos no son lineales, así que restar una tendencia simplemente lineal puede modificar la señal de manera apreciable según el caso.
Acerca del tema del residuo, la descomposición ejemplificada en la figura19, que correspon- de al uso del método EEMD, parece peculiar si se compara con las demás. El residuo tiene un aspecto muy similar a las dos IMFs anteriores a él, y, en comparación con los residuos obtenidos con el método EMD y con el método CEEMDAN, su valor es mucho menor. Un sencillo procedimiento que se puede aplicar para analizar la calidad de la descomposición con el método EEMD es el que se muestra en la figura29. En las tres figuras se representa la serie temporal estudiada como ejemplo hasta ahora (14, izquierda), a la que se le ha resta- do el residuo de cada descomposición, de manera similar a lo que hacen Kuridze et al. [4].
También se representa la IMF de cada método que caracteriza la oscilación de la espícula. Es evidente que el método EEMD ofrece un residuo que no corresponde a la tendencia real de los datos. Los métodos EMD y CEEMDAN presentan una concordancia mucho mayor entre la señal sin tendencia y la IMF característica de la oscilación. Comparando estos dos méto- dos, el método CEEMDAN se destaca, ya que la IMF es muy similar a la señal sin tendencia.
Cabe recordar que el método EMD descompone la señal en más IMFs en este caso, y pre- senta mezcla de modos. Todo esto puede contribuir a la mayor discrepancia que presenta en comparación con el método CEEMDAN.
Figura 29: Comparación entre restar el residuo a la serie temporal del panel izquierdo de la figura14(azul) con la IMF que refleja la oscilación de la espícula (naranja). Arriba: Método EEMD (izquierda) y método CEEMDAN (derecha). Abajo: Método EMD.
.
Debido a la manera en la cual el método EEMD obtiene las IMFs, con un promedio sobre un colectivo deN realizaciones de la misma descomposición, es posible que en diferentes rea- lizaciones se obtenga un número de IMFs distinto, y esto puede afectar al resultado cuando se hace el promedio. Esta es una posible explicación a la gran similitud que hay entre las dos últimas IMFs y el residuo de la descomposición presentada en la figura19.
5. Conclusiones y Observaciones
En este trabajo se ha presentado un nuevo método para la descomposición de funciones pe- riódicas y dos de sus variantes. El análisis ‘clásico’ de Fourier tiene graves limitaciones desde un punto de vista científico, ya que impone condiciones como la linealidad de los sistemas.
En ciencia, en particular en Física, pocos sistemas reales son lineales o estacionarios. En lu- gar de utilizar el análisis de Fourier, aquí se propone el uso de la Descomposición Empírica de Modos, que no impone ninguna de estas condiciones.
En la EMD, se descompone una señal en Funciones Intrínsecas de Modo (IMF), que son mu- cho más representativas de los diferentes modos de oscilación presentes en una señal que las funciones trigonométricas utilizadas en el análisis de Fourier. También se obtiene una descomposición más simple, ya que el número de autofunciones es reducido en compara- ción con el gran número de términos necesarios para reconstruir adecuadamente la señal en función de senos y cosenos. El método EMD, además, es guiado por los datos; las IMFs son autofunciones que cambian según la señal que se descomponga, no vienen prefijadas.
El mayor inconveniente del método EMD deriva de la localidad de la definición de las IMFs, algo que resulta en la mezcla de modos que se ha observado. Debido a esto, mejoras en el algoritmo para la obtención de IMFs se han presentado y estudiado. El método EEMD in- troduce ruido blanco aleatorio a la señal, resultando en un colectivo de descomposiciones que luego se promedian. Sin embargo, no soluciona el problema de la mezcla de modos y además la descomposición ‘promedio’ presenta un número inusual de IMFs en compara- ción a los otros métodos. La IMF característica de la oscilación que se quiere estudiar puede verse comprometida por el ruido, como se observa en la figura21al compararla con la IMF correspondiente al método EMD original. La amplitud de oscilación se ve afectada por la nueva descomposición, lo que puede llevar a errores en el análisis posterior.
También se introduce el método CEEMDAN, que utiliza también ruido para solucionar el problema de la mezcla de modos, pero de manera distinta al método EEMD. Esto resulta no solo en la solución al problema de la mezcla de modos en las zonas alejadas de la frontera, donde hay problemas con la transformada de Hilbert y la mezcla es difícil de evitar, sino que el problema se soluciona sin afectar a los datos. Al comparar las IMFs características para una misma señal temporal con los tres métodos distintos (figura25), el método CEEMDAN no parece afectar a la amplitud de la misma manera que el método EEMD.
Insistiendo en el tema del residuo en la comparación de los tres algoritmos estudiados para la obtención de las IMFs, el residuo obtenido con los algoritmos EMD y CEEMDAN parece reflejar mucho más la tendencia que presenta la espícula que no el que resulta de aplicar el método EEMD.
Se han reproducido los cálculos hechos por Kuridze et al. [4], donde utilizan la descompo- sición en senos y cosenos (análisis de Fourier) para obtener un período de oscilación pro- medio de 180 segundos. Sin embargo, estos autores no solo restan una tendencia lineal a los datos (los residuos de las descomposiciones en IMFs no son lineales), sino que como se ha observado la oscilación de la espícula no se asemeja a la de una función sinusoidal simple, es mucho más compleja puesto que ni su amplitud ni su frecuencia instantáneas son cons- tantes. Con la descomposición en IMFs se obtiene una función que refleja bastante mejor la naturaleza real de la oscilación transversal de la espícula solar. El período de oscilación que ellos obtienen puede caracterizar adecuadamente la oscilación de la espícula (en las descomposiciones con EMD los valores del período son cercanos a 200 segundos), pero el análisis con EMD refleja mucho mejor el carácter no estacionario de la oscilación.
Referencias
[1] Marcelo Alejandro Colominas. “Métodos guiados por los datos para el análisis de seña- les: contribuciones a la descomposición empírica en modos”. Tesis doct. Universidad Nacional del Litoral, 2016.
[2] Herman H. Goldstine.A History of Numerical Analysis From the 16th Through the 19th Century. Springer, 1976.
[3] NE Huang y col. “The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis”. En:Proceedings of the Royal So- ciety A: Mathematical, Physical and Engineering Sciences(1998). ISSN: 13645021. DOI: 10.1098/rspa.1998.0193.
[4] D. Kuridze y col. “Characteristics of transverse waves in chromospheric mottles”. En:
Astrophysical Journal(2013).ISSN: 15384357.DOI:10.1088/0004-637X/779/1/82. arXiv:
1310.3628[astro-ph.SR].
[5] Juan Martínez-Sykora y col. “Two-dimensional Radiative Magnetohydrodynamic Simu- lations of Partial Ionization in the Chromosphere. II. Dynamics and Energetics of the Low Solar Atmosphere”. En:The Astrophysical Journal847 (2017).ISSN: 1538-4357.DOI: 10.3847/1538-4357/aa8866. arXiv:1708.06781[astro-ph.SR].
[6] Solar y Astrophysics Lab.Scientists Explain Mysterious Plasma Jets On The Sun. 2004.
URL:http://www.lmsal.com/Press/spicules2004/.
[7] María E Torres y col. “A Complete Ensemble Empirical Mode Decomposition with Adap- tive Noise”. En:IEEE International Conference on Acoustics, Speech and Signal Proces- sing (ICASSP)(2011).ISSN: 15206149.DOI:10.1109/ICASSP.2011.5947265. arXiv:arXiv:
1011.1669v3.
[8] ZHAOHUA WU y NORDEN E. HUANG. “ENSEMBLE EMPIRICAL MODE DECOMPOSI- TION: A NOISE-ASSISTED DATA ANALYSIS METHOD”. En:Advances in Adaptive Data Analysis(2009).ISSN: 1793-5369.DOI:10.1142/S1793536909000047.
