Cette communication porte sur l'estimation des paramètres de modes vibratoires ayant des excitations corrélées. L'application visée est l'analyse des modes vibratoires du soleil. Deux algorithmes d'estimation sont proposés. Le premier, qui est une extension de la méthode habituellement utilisée en héliosismologie, repose sur les propriétées asymptotiques du périodogramme. Cette méthode présente l'avantage de permettre une prise en compte immédiate des contraintes sur le profil des modes. Le second repose sur une modélisation temporelle du signal par une somme bruitée de processus autorégressifs ayant des excitations corrélées. Après une estimation des paramètres dynamiques la matrice de covariance des excitations est estimée par un algorithme EM. Les performances des deux méthodes sont évaluées par des simulations de Monte Carlo et comparées à la borne de Rao Cramèr. Enfin, des résultats sur des signaux mesurés par le satellite SOHO sont présentées

CERCOS, Miguel

FERRARI, André

GELLY, Bernard

SCHMIDER, Francois Xavier

I-Revues

Estimation Spectrale Paramétrique de Modes Corrélés,
Application à l'Héliosismologie
André Ferrari, François Xavier Schmider, Miguel Cercos, Bernard Gelly
UMR 6525 Astrophysique
Université de Nice Sophia-Antipolis, Parc Valrose,
06108 Nice CEDEX 2, FRANCE
ferrari@unice.fr
Résumé  Cette communication porte sur l'estimation des paramètres de modes vibratoires ayant des excitations corrélées.
L'application visée est l'analyse des modes vibratoires du soleil. Deux algorithmes d'estimation sont proposés. Le premier, qui
est une extension de la méthode habituellement utilisée en héliosismologie, repose sur les propriétées asymptotiques du pério-
dogramme. Cette méthode présente l'avantage de permettre une prise en compte immédiate des contraintes sur le prol des
modes. Le second repose sur une modélisation temporelle du signal par une somme bruitée de processus autorégressifs ayant des
excitations corrélées. Après une estimation des paramètres dynamiques la matrice de covariance des excitations est estimée par
un algorithme EM. Les performances des deux méthodes sont évaluées par des simulations de Monte Carlo et comparées à la
borne de Rao Cramèr. Enn, des résultats sur des signaux mesurés par le satellite SOHO sont présentées.
Abstract  The purpose of this communication is the parameter estimation of vibration modes having correlated excitations.
The focused application are the solar vibration modes. Two estimation algorithm are proposed. The rst in an extension of the
method classically applied in helisismology and relies on the asymptotic properties of the periodogram. The advantage of this
method is that it easily takes into account constraints on modes proles. The second relies on the temporal modeling of the signal
has a noisy sum of autoregressive processes having correlated excitations. After the estimation of the dynamical parameters, the
excitations covariance matrix is estimated by an EM algorithm. The performances of the methods are evaluated by Monte Carlo
simulations and compared to the Cramèr Rao bound. Finally, results obtained from SOHO satellite signal are presented.
1 Introduction
Cette communication porte sur l'analyse spectrale de
données héliosismiques : l'onde stationnaire piégée dans
une cavité à l'intérieur du soleil induit des variations de la
vitesse radiale à la surface. L'amplitude de ces variations
vérie ainsi une équation diérentielle stochastique du se-
cond ordre. La gure 1 représente le périodogramme du
signal mesuré par l'instrument GOLF du satellite SOHO.
L'objectif est d'analyser ce signal dans des bandes de fré-
quences contenant des modes de même degré.
L'analyse de ce signal a été jusqu'à présent réalisée en
supposant que les excitations des diérents modes étaient
indépendantes. La densité spectrale de puissance (DSP)
est alors modélisée par une somme de Lorentziennes dont
les paramètres sont estimés par maximum de vraisem-
blance à partir du périodogramme, [4].
Nous proposons dans cette communication de lever l'hy-
pothèse de décorrélation des excitations. Pour cela, deux
méthodes d'estimation prenant en compte leurs coecients
de corrélation sont proposées. La première est une exten-
sion de la technique actuellement utilisée en héliosismo-
logie. La seconde, qui constitue la contribution majeure
de cette communication, est basée sur la représentation
de chaque mode par un processus autorégressif du se-
cond ordre. Ces deux méthodes seront référencées dans la
suite par approche fréquentielle et temporelle. Après
avoir présenté les deux algorithmes, leurs performances
sont étudiées par des simulations et comparées à la borne
de Rao Cramèr.
2 Méthodes d'estimation
2.1 Modèle spectral
Chaque mode est individuellement solution d'une équa-
tion diérentielle stochastique du second ordre. On im-
posera la contrainte, habituelle en héliosismologie, que à
l'intérieur d'une bande de fréquence la puissance de l'ex-
citation ainsi que les amortissements des diérents modes
sont constants.
Sous ces hypothèses et compte tenu de l'important fac-
teur de sur-échantillonnage, la DSP du signal bruité s'écrit :
S(f) = 
2
e
(
L
X
k=1
jH
k
(f)j
2
+
L
X
k;l=1
r
kl
H
k
(f)H
l
(f)

) + 
2
w
;
(1)
avec :
H
k
(f) = 1=(1  (f=f
k
)
2
+ j2(f=f
k
)); (2)
où 
2
e
désigne la puissance des excitations et r
kl
, (jr
kl
j < 1)
leurs coecients de corrélation. La gure 2 montre l'eet
d'une variation de r
12
sur S(f) dans le cas L = 2 pour :
f
1
= 9:2, f
2
= 10:8,  = 0:015, 
2
e
= 0:01 and 
2
w
= 0:5
Le processus mesuré étant supposé linéaire et la réso-
lution fréquentielle recherchée nécessitant un très grand
nombre d'échantillons ( 10
6
), le périodogramme dans la
0 1 2 3 4 5 6 7
0
2
4
6
8
10
12
x 104
D
en
si
te
 s
pe
ct
ra
le
Frequence (mHz)
Fig. 1: Périodogramme de 600 jours, T
e
= 20s. Chaque
pic est un multiplet de 1,2 ou 3 composantes.
6 7 8 9 10 11 12 13 14
0
2
4
6
8
10
12
14
S(
f)
f
r12=0.9
r12=0.5
r12=0
r12=−0.5
r12=−0.9
Fig. 2: Eet du coecient de corrélation sur la DSP pour
un mode double.
bande à analyser va vérier les propriétés asymptotiques
classiques d'indépendance et de distribution suivant une
loi S(f)
2
2
, [6]. On propose donc d'estimer le vecteur de
paramètres :
 = f
2
e
;
2
w
; ; (f
k
)
k=1::: ;L
; (r
kl
)
k=1:::L;
l=1;::: ;k 1
g (3)
par maximum de vraisemblance :
^
 = argmin

N
f
X
k=N
i
(lnS(f
k
) + I(f
k
)=S(f
k
));
où les I(f
k
), k = N
i
; : : : ; N
f
désignent les valeurs du pé-
riodogramme des mesures sur l'intervalle de fréquence à
analyser.
2.2 Modèle temporel
Nous allons supposer dans cette partie que chaque mode
peut être représenté par un processus autorégressif du se-
cond ordre. Pour L modes, les mesures y
n
peuvent donc
s'écrire :
x
(k)
n
=  a
(k)
1
x
(k)
n 1
  a
(k)
2
x
(k)
n 2
+ e
(k)
n
; (4)
y
n
=
L
X
k=1
x
(k)
n
+ w
n
; n = 0; : : : ; N   1; (5)
où e
(k)
n
et w
n
sont des bruits blancs. Il est important de
souligner que les contraintes sur l'égalité des amortisse-
ments imposées dans le modèle précédent ne peuvent être
utilisées ici, le signal étant obtenu en ramenant en bande
de base la plage de fréquence étudiée.
L'ensemble des mesures sera noté Y
N
= fy
0
; : : : ; y
N
g.
On notera,
~
E
n
= (e
1
n
; : : : e
L
n
)
t
et Q = E [
~
E
n
~
E
t
n
]. Le coef-
cient de corrélation de l'excitation des modes i et j est
alors r
ij
= Q
ij
=
p
Q
ii
Q
jj
. Le problème est d'estimer les
coecients autoregréssifs notés , la matrice de covariance
Q et la puissance du bruit de mesure 
2
w
.
La première partie de l'algorithme consiste à estimer .
Considérons pour cela le signal z
n
obtenu en ltrant y
n
par:
H(z) =
L
Y
q=1
(1 + a
(q)
1
z
 1
+ a
(q)
2
z
 2
): (6)
Il est aisé de vérier que E z
n
z
n+l
= 0 pour jlj > 2L. z
n
possède donc une représentation MA et en conséquence,
y
n
une représentation ARMA où les coecients du déno-
minateur sont ceux de H(z).
On propose donc d'estimer H(z) en appliquant à y
n
un algorithme d'estimation ARMA. On utilise pour cela
l'algorithme d'estimation par sous-espaces, [8, p.115]. Les
ltres 1 + a
(k)
1
z
 1
+ a
(k)
2
z
 2
son ensuite obtenus en com-
binant les L paires de racines conjuguées de H(z).
La seconde partie de l'algorithme consiste à remplacer
 par son estimée
^
 an d'estimer Q et 
2
w
. Une solution
ecace de ce problème nécessite un modèle faisant inter-
venir de façon explicite le vecteur
~
E
n
. Nous proposons
pour cela d'utiliser un modèle d'état. Chaque mode peut
se modéliser par :
~
X
(k)
n+1
= F
(k)
~
X
n
+G
(k)
~
E
(k)
n
; x
(k)
n
= H
(k)
t
~
X
n
; (7)
avec :
F
(k)
=
 
0  a
(k)
2
1  a
(k)
1
!
; G
(k)
= H
(k)
=

0
1

: (8)
Ainsi, un modèle d'état de y
n
est :
~
X
n+1
= F
~
X
n
+G
~
E
n
; y
n
= H
t
~
X
n
+ w
n
; (9)
où F est une matrice 2L  2L bloc diagonale avec les
F
(k)
sur la diagonale. H est un vecteur 2L  1 obtenu
par concatenation des H
(k)
et G est une matrice 2L  L
avec pour éléments non nuls G(2k; 2) = 1, k = 1; : : : ; L.
L'ensemble des états sera noté X
N
= f
~
X
0
; : : : ;
~
X
N
g.
Nous allons considérer les estimations de Q et 
2
w
ob-
tenues en maximisant L(Y
N
;
^
;Q; 
2
w
) où L(Y
N
; ;Q; 
2
w
)
est la log-vraisemblace des observations. Nous utiliserons
pour cela l'algorithme EM. Dans ce cas, Y
N
va jouer le rôle
de données incomplètes et X
N
de données complètes.
Cet algorithme a été utilisé pour l'estimation de l'ensemble
des paramètres d'un modèle d'état dans [1]. Nous présen-
terons ici pour le cas où seule la covariance de l'excitation
et la puissance du bruit d'observation sont estimées. On
pourra se référer à [3] pour plus de détails.
L'étape E réside dans le calcul de l'espérance de la log-
vraisemblance de X
N
conditionnée à Y
N
et les valeurs cou-
rantes de Q et 
2
w
, notées
^
Q
[k]
, 
2;[k]
w
.
0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
0.08
0.09
0.1
0.11
0.12
0.13
0.14
0.15
∆ f
BR
C(
r 12
)
Fig. 3: BRC sur r
12
en fonction de f = f
2
  f
1
pour un
mode double.
6 7 8 9 10 11 12 13 14
0
2
4
6
8
10
12
f
D
SP
 m
od
es
Fig. 4: DSP correspondant aux deux valeurs extrêmes de
f pour le calcul de la BCR.
Dans notre cas, la log-vraisemblace de X
N
se calcule fa-
cilement en utilisant la propriété Markovienne du modèle
d'état:
L(X
N
;Q) =  
N   1
2
(log jQj+ log
2
w
) 
1
2
N 1
X
k=1
(G
t
~
V
t
k
Q
 1
~
V
k
G+ w
2
k
=
2
w
);
~
V
k
=
~
X
k
  F
~
X
k 1
; w
k
= y
k
 H
t
~
X
k
:
La maximisation de l'espérance conditionelle de l'ex-
pression précédente (étape M) donne alors :
^
Q
[k+1]
=
1
N
N 1
X
k=1
G
t
E [
~
V
k
~
V
t
k
=Y
N
;
^
Q
[k]
; 
2;[k]
w
]G; (10)

2;[k+1]
w
=
1
N
N 1
X
k=1
E [w
2
k
=Y
N
;
^
Q
[k]
; 
2;[k]
w
]: (11)
Les quantités intervenant dans la somme peuvent être cal-
culées en utilisant un algorithme de lissage à interval xe,
[1]. On peut souligner que la structure de (10) garantie que
^
Q
[k]
est dénie semi-positive. Par ailleurs, les estimateurs
par maximum a posteriori de Q est 
2
w
sont, [7, 2]:
^
Q
MAP
=
1
N   1
N 1
X
k=1
G
t
^
~
V
k=N
^
~
V
t
k=N
G; (12)

2
w;MAP
=
1
N
N 1
X
k=1
w
2
k
: (13)
Une implantation itérative de ces estimateurs s'obtient
en calculant les vecteurs d'état lissés à l'itération k à par-
tir de
^
Q
[k]
MAP
et 
2;[k]
w;MAP
et en calculant
^
Q
[k+1]
MAP
et 
2;[k+1]
w;MAP
à partir de (12,13). Cet algorithme correspond à l'algo-
rithme EM décrit plus haut si l'on considère que la mise à
jour de Q et 
2
w
se fait en supposant que les variance de
^
~
X
k=n
et les covariance entre
^
~
X
k=n
et
^
~
X
k 1=n
sont nulles.
Un défaut majeur de cette solution par rapport à EM est
bien sur qu'elle ne garantie pas une augmentation de la
vraisemblance à chaque itération.
2.3 Bornes de Rao Cramèr
Le but de cette section est d'étudier les bornes de Rao
Cramèr (BRC) des paramètres caractéristiques du modèle
de modes corrélés contenu dans une bande de fréquence
[f
m
; f
M
]. Nous proposons d'utiliser pour celà l'expression
asymptotique (N grand) des termes de la matrice d'in-
formation de Fisher d'un processus Gaussien en fonction
de sa DSP, [5, p. 77]. Si l'on modélise la DSP par (1), le
terme I
k;l
de la FIM correspondant aux composantes k et
l de (3) est alors donné par :
I
k;l
=
N
2(f
M
  f
m
)
Z
f
M
f
m
1
S
2
(f)
@S(f)
@
i
@S(f)
@
j
df: (14)
Une fois les dérivées partielles calculées, les termes de la
FIM sont alors obtenus par intégration numérique en uti-
lisant la fonction quad8.m de MATLAB.
La gure 3 donne la BRC sur r
12
pour f
1
= 10 f=2,
f
2
= 10+f=2, r
12
= 0:5,  = 0:015, 
2
e
= 0:01, 
2
w
= 0:5
et N = 1024 en fonction de f . La gure 4 représente
les DSP correspondant aux deux valeurs extrêmes de f .
L'allure concave de cette courbe est assez intuitive. Si f
est petit, les deux excitations sont mélangées dans la même
bande de fréquence et il est donc dicile de les séparer
et si f est trop grand, les excitations sont décorélées
par les ltres. Cette courbe montre également les faibles
performances auxquelles on peut s'attendre. En eet le
plus faible écart type que l'on peut obtenir pour r
12
dans
ce cas est 0.3 : 15% environ de l'intervalle de variation de
r
12
.
3 Simulations
Cette première simulation a pour objectif l'évaluation
des performances des approches temporelles et fréquen-
tielles par une méthode de Monte Carlo dans le cas L = 2.
An de pouvoir comparer les résultats, les paramètres des
deux modèles ont été choisis de façon à ce que les DSP
de chacun des modes aient dans les deux cas même : fré-
quence maximale, amplitude maximale et largeur à 3dB.
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
0
20
40
60
80
100
120
f
D
SP
 M
od
èl
es
 A
R
50 50.05 50.1 50.15 50.2 50.25 50.3 50.35 50.4 50.45 50.5
0
20
40
60
80
100
120
f
D
SP
 M
od
èl
es
 fr
éq
ue
nt
ie
lle
s
Fig. 5: DSP des modes pour les modèles autorégressifs et
fréquentiels.
Tab. 1: Performance des deux approches.
r
12
(r
12
) f
1
(f
1
) f
2
(f
2
)
temp. -0.46 0.14 0.20 3.2e-3 0.30 3.2e-3
fréq. -0.53 0.16 0.20 2.1e-3 0.30 2.2e-3
BRC - 0.18 - 2.0e-3 - 2.0e-3
Dans le second cas, f
m
= 50 et f
M
= 50; 5. La gure 5
présente les deux modes correspondants aux deux modèles
du signal. Le tableau 1 présente les moyennes, écarts types
et racines carrées de la BRC pour r
12
, f
1
et f
2
obtenues à
partir de 256 réalisations d'un signal de 500 échantillons.
Pour l'estimation des fréquences, les performances de la
méthode fréquentielle sont légérement supérieures à celles
de la méthode temporelle. Les résultats sont en accord
avec la BRC. Cette tendance est inversée pour l'estimation
du coecient de corrélation. Il est important de constater
que les écarts types obtenus sont supérieurs à ceux donnés
par la BRC, ce qui s'explique par le fait que, comme le
montre la première colonne du tableau, les estimées de
r
12
sont légèrement biaisées dans les deux méthodes.
Enn, les deux méthodes ont été appliquées à un signal
de 600 jours échantillonné à 20s provenant de l'instrument
SOHO. Pour la méthode temporelle, le signal correspon-
dant aux bandes de fréquence analysées a été démodulé et
décimé par transformée de Fourier directe et inverse. Le
tableau 2 présente les résultats obtenus pour diérentes
bandes. Le nombre d'échantillons correspondants à chaque
bande est environ 700. On constate, contrairement aux fré-
quences, une très mauvaise adéquation des résultats obte-
nus pour l'estimation du coecient de corrélation.
La gure 6 présente les deux DSP estimées ainsi que le
périodogramme pour la bande [1746; 1752]Hz. Un fac-
teur essentiel dans le résultat obtenu est le fort a priori
imposé sur l'amplitude des modes dans l'approche fré-
quentielle. Cette contrainte qui a peu d'eet sur l'esti-
mation des fréquences a une inuence fondamentale sur le
coecient de corrélation. Le problème crucial est donc la
validation de cette hypothèse qui au vu de la gure 6 ne
semble pas acquise.
Tab. 2: Résultats obtenus sur le signal Golf.
bande (Hz) r^
12
^
f
1
^
f
2
[1611,1614] temp. -0.38 1612.32 1613.22
fréq. 0.96 1612.32 1613.14
[1746,1752] temp. -0.63 1748.85 1749.79
fréq. -0.65 1748.85 1749.70
[1882,1888] temp. -0.75 1884.76 1885.55
fréq. 0.44 1884.71 1885.48
1746 1747 1748 1749 1750 1751 1752
0
5
10
15
20
25
30
35
f
D
SP
periodogramme
frequentielle
temporelle
Fig. 6: Résultat sur un mode double du signal GOLF.
Références
[1] L. Deng and X. Shen. Maximum likelihood in statis-
tical estimation of dynamic systems: Decomposition
algorithm and simulation results. EURASIP Signal
Processing, 57:6579, 1997.
[2] G. Favier, G. Alengrin, and P. Guillermin. Algorithmes
de ltrage adaptatif, estimation des covariances de
bruit dans un système linéaire stochastique. Technical
report, Laboratoire de Signaux et Systèmes, 1979.
[3] A. Ferrari, F.-X. Schmider, B. Gelly, and G. Alengrin.
Parametric representation of helioseismic data. part 1:
a method for estimation of mode correlation. Astro-
nomy and Astrophysics, 1999. à paraître.
[4] D. Fierry Fraillon. Étude des densités spectrales
des modes acoustiques solaires: une modélisation des
modes de bas degrés. PhD thesis, Université de Nice
Sophia-Antipolis, 1999.
[5] S. M. Kay. Fundamentals of Statistical Signal Pro-
cessing : Estimation Theory. Signal Processing Series.
Prentice Hall, 1993.
[6] M. B. Priestley. Spectral Analysis And Time Series.
Monographs and Textbooks. Academic Press, 8th edi-
tion, 1994.
[7] A.P. Sage and J.L. Melsa. System Identication, vo-
lume 80. Academic Press, 1971.
[8] Petre Stoica and Randolph Moses. Introduction to
Spectral Analysis. Prentice Hall, 1997.


Estimation spectrale paramétrique de modes corrélés, application à l'héliosismologie

Abstract

Similar works

Full text

Available Versions

I-Revues