FR2893434A1 - Procede et dispositif d'identification de parametre multidimensionnels : application a la localisation et la reconstruction d'activites electriques de profondeur au moyen d'observations de surfaces - Google Patents

Abstract

Procédé d'identification de paramètres multidimensionnels propres à une pluralité de P>=1 sources d'intérêt présentes au sein d'un environnement conducteur multidimensionnel prédéterminé, au moyen d'une pluralité d'observations (60) en nombre N>=1 fini.Selon l'invention, on exploite i) la factorisation de la formulation matricielle du problème, ii) la création d'un réseau virtuel d'ordre 2q (q>1) de capteurs de par l'exploitation de cumulants d'ordre 2q des observations, et iii) le concept de déflation étendu à l'ordre 2q en tenant compte de la présence de sources potentiellement (mais pas totalement) corrélées.Application du procédé et du dispositif selon l'invention aux domaines de :- l'électroencéphalographie ;- la magnétoencéphalographie ;- la géophysique ;- la sismologie.

Description

11110 r. 1

Procédé et dispositif d'identification de paramètres multidimensionnels : application à la localisation et la reconstruction d'activités électriques de profondeur au moyen d'observations de surface 1. Domaine de l'invention Le domaine de l'invention est celui de l'acquisition et du traitement de signaux représentatifs d'activités générées par un ensemble de sources internes à un environnement multidimensionnel donné que l'on cherche à étudier. Plus précisément, l'invention concerne une nouvelle technique d'identification de paramètres multidimensionnels qui permet, entre autres, la localisation et la reconstruction d'activités électriques, appelées communément sources, générées au sein d'un environnement multidimensionnel, à partir des seules observations acquises en certains points dudit environnement au moyen d'un ensemble de capteurs physiques. L'invention s'inscrit donc dans un contexte et une problématique technique pouvant intéresser de nombreuses disciplines pour lesquelles ce problème d'identification de paramètres multidimensionnels de sources d'intérêt à partir d'observations est primordial et nécessaire à la meilleure compréhension des phénomènes internes à l'environnement étudié. Ces disciplines concernent à titre d'exemple illustratif et non limitatif, aussi bien le domaine biomédical et plus précisément l'électrophysiologie humaine ou animale, que le domaine de la géophysique, par exemple pour l'estimation de la position de l'épicentre et de la propagation des séismes. 2. État de la technique 2.1 Contexte spécifique à l'électrophysiologie Dans les domaines complémentaires de l'ElectroEncéphaloGraphie (EEG) et de la MagnétoEncéphaloGraphie (MEG), les activités post-synaptiques des neurones corticaux peuvent être enregistrées de façon totalement non invasive, c'est-à-dire sans besoin d'implanter des électrodes intracérébrales. L'EEG recueille l'activité électrique de surface au moyen de capteurs 30 (électrodes d'enregistrement) disposés de façon standardisée sur le cuir chevelu et permet de représenter l'évolution dans le temps d'une différence de potentiels électriques entre chaque électrode et une référence commune. Parallèlement, la MEG enregistre le champ magnétique (d'amplitude très faible) induit par l'activité électrique cérébrale grâce à des capteurs ultrasensibles 5 disposés sur un casque couvrant la totalité du scalp. Chez le sujet sain, ces deux techniques permettent d'enregistrer des rythmes spontanés particuliers liés à des états physiologiques normaux (veille, sommeil, attention) ; elles permettent également de recueillir les réponses évoquées par une stimulation donnée, qui reflètent l'activation des différentes 10 structures cérébrales mises en jeu lors de processus perceptuels et/ou cognitifs. De façon similaire, les enregistrements EEG et MEG permettent de mieux comprendre les mécanismes impliqués dans le cadre de certaines pathologies neurologiques telles que l'épilepsie, la maladie d'Alzheimer ou les tumeurs cérébrales. 15 Dans tous les cas, l'analyse des bases neuronales d'un système cognitif complexe normal ou d'un dysfonctionnement pathologique, tel qu'une crise d'épilepsie par exemple, nécessite l'identification précise à la fois des régions cérébrales participantes et de la séquence temporelle d'activation entre ces régions. 20 Les techniques connues d'imagerie fonctionnelle (IRMf pour Imagerie par résonance magnétique fonctionnelle, TEP pour Tomographie à Emission de Positons, TEMP pour Tomographie par Emission Mono Photonique) - qui bénéficient le plus souvent d'une haute résolution spatiale (de l'ordre du mm) - ont prouvé leur efficacité dans la définition des aires anatomiques activées lors 25 d'opérations cognitives. En revanche ces méthodes restent sévèrement limitées dans leur capacité à révéler l'information temporelle de cette activation puisqu'au mieux, elles fournissent une image d'activation moyenne sur plusieurs centaines de millisecondes. Inversement, l'EEG et la MEG, de par leur excellente résolution 30 temporelle 1 ms), autorisent une étude fine des dynamiques d'activation neuronale. Cependant, l'identification précise des régions cérébrales activées nécessite l'utilisation d'outils mathématiques à la fois fiables et précis permettant de résoudre le problème inverse, c'est-à-dire de localiser et de reconstruire les sources de l'activité cérébrale considérée à partir des seules observations de surface.

Comme illustré sur la figure 1, le problème inverse en électrophysiologie humaine se définit comme la possibilité, à partir des simples enregistrements de surface (potentiels électriques et/ou champs magnétiques) 11 de l'activité cérébrale issus de capteurs 10 et en utilisant des modèles de tête 12 et de sources appropriés au milieu de conduction considéré, d'identifier les régions cérébrales 13 responsables des activités EEG et/ou MEG enregistrées. Plus particulièrement, le problème inverse en EEG et en MEG consiste à estimer les paramètres des sources dipolaires de l'activité cérébrale et à reconstruire les décours temporels associés, à partir des seules observations de surface. De manière générale, la résolution du problème inverse en EEG et en MEG nécessite que l'on établisse : un modèle de sources, qui doit rendre compte des caractéristiques spatiales et temporelles des sources neuronales à l'origine de l'activité électromagnétique cérébrale ; et un modèle de volume conducteur, qui doit reproduire au mieux la géométrie et les propriétés physiques de tous les constituants de la tête ; et que l'on résolve le problème direct associé, qui vise à caractériser la conduction de l'activité des sources au sein du volume conducteur. Ces différents points sont décrits plus précisément ci-dessous. Notons que le problème inverse n'est pas spécifique à l'électrophysiologie, puisqu'on le rencontre également, par exemple, dans le domaine de la sismologie, où l'on cherche à estimer l'épicentre et la propagation des séismes, tel que décrit dans l'ouvrage de A. TARANTOLA : Inverse Problem Theory et Model Parameter Estimation , réédité par SIAM en 2004. 2.1.1 Modèle de sources Pour représenter l'activité électromagnétique cérébrale, le modèle de source le plus couramment utilisé est le modèle dipolaire, qui assimile l'activité d'une petite zone corticale à celle d'un dipôle de courant, tel que représenté sur la figure 2. On distingue alors deux cas de figure. Dans le premier cas, le nombre de dipôles est supposé être inférieur ou égal au nombre d'observations de surface (activités EEG ou MEG). On parle alors de mélange surdéterminé de sources, ou de problème bien posé . Cependant, il est en pratique inexact de réduire le fonctionnement cérébral, normal ou pathologique, à l'activité d'un nombre restreint de sources. On considère alors un second cas de figure, dans lequel le nombre de sources est supposé être strictement supérieur au nombre d'observations de surface. Le mélange de sources est alors dit sous-déterminé, et le problème est mal posé . Dans ce second cas de figure, il est alors important de distinguer le problème de la localisation des dipôles de celui de la reconstruction des activités électriques de profondeur générées par ces derniers. En effet, alors que le second problème ne peut théoriquement pas être résolu de manière unique sans l'ajout et l'exploitation d'information a priori sur les sources d'intérêt, il en est tout autrement pour le premier problème. Les inventeurs ont constaté que ce résultat est méconnu des chercheurs du domaine biomédical. 2.1.2 Modèle de volume conducteur Outre les caractéristiques propres des sources, le potentiel électrique et le champ magnétique enregistrés à la surface du scalp d'un sujet sont également liés aux contraintes physiques et géométriques des différents tissus de la tête. La tête est assimilée à un volume conducteur qui doit tenir compte des inhomogénéités des différents milieux (cerveau, liquide céphalo-rachidien, crâne et peau). La tête est donc généralement modélisée par un ensemble de trois ou quatre couches concentriques, de conductivités différentes, représentant les différents tissus traversés lorsque le signal issu des sources atteint la surface cutanée.

Les conductivités de chacun de ces milieux peuvent être considérées comme isotropes, tel que décrit dans l'article de S. RUSH et D. A. DRISCOLL, EEG electrode sensitivity ù An application of reciprocity , IEEE Transactions On Biomedical Engineering, vol. 38, pp. 15ù22, Janvier 1969, ou bien encore anisotropes, comme décrit dans l'article de J. C. De MUNCK et M. J. PETERS, A fast method to compute the potential in the multisphere model , IEEE 5 Transactions On Biomedical Engineering, vol. 40, pp. 1166ù1174, Novembre 1993. Le modèle de géométrie le plus simple est le modèle sphérique (figure 3.a) qui assimile la tête à un ensemble de sphères concentriques, où chaque couche correspond à un tissu différent. Le modèle le plus couramment utilisé est un modèle à trois sphères qui représentent respectivement le cerveau 30, le crâne 31 et le scalp 32 d'un sujet. Le modèle sphérique n'est cependant qu'une approximation grossière de la géométrie de la tête. Des modèles à géométrie plus réaliste (figure 3.b) ont donc été développés, et sont construits, pour chaque sujet, à partir d'images IRM anatomique. Des méthodes permettent en effet d'extraire, par segmentation des images IRM, les contours des trois structures d'intérêt que sont le cerveau 33, le crâne 34 et le scalp 35, puis de générer des maillages 3D de ces trois surfaces. 2.1.3 Résolution du problème direct Résoudre le problème direct en EEG et en MEG consiste à savoir calculer 20 le champ électromagnétique généré à la surface du scalp par une configuration de sources connues dans le volume cérébral. L'application des lois de la physique telles que les équations de Maxwell, la loi de conservation des charges et la loi de Biot et Savart, permet de calculer le potentiel électrique et le champ magnétique créés sur les capteurs de surface, 25 connaissant la configuration des sources intracérébrales, et la géométrie et les conductivités des différents tissus de la tête. Le choix de la méthode de calcul du potentiel électrique et du champ magnétique va dépendre du type de modèle de tête. Dans le cas d'un modèle sphérique (figure 3.a), on peut calculer analytiquement le potentiel électrique et le 30 champ magnétique engendrés par un dipôle donné, tel que décrit respectivement dans les articles de P. BERG et SHERG, A fast method for forward computation of multiple-shell spherical head models , Electroencephalography and Clinical Neurophysiology, vol. 90, no. 1, pp. 58-64, Janvier 1994, et de J. SARVAS, Basic mathematical and electromagnetic concepts of the biomagnetic inverse problems , Physics in Medicine and Biology, vol. 32, pp. 11-22, 1987.

Dans le cas des modèles réalistes précités (figure 3.b), le calcul du potentiel électrique et du champ magnétique n'est réalisable qu'à l'aide de méthodes numériques connues de l'art antérieur, à savoir : la méthode des éléments frontières (BEM pour Boundary Element Method, en anglais), bien décrite dans l'article de A. S. FERGUSON, X. ZHANG, et G. STROINK, "A complete linear discretization for calculating the magnetic field using the boundary element method," IEEE Transactions On Biomedical Engineering, vol. 41, pp. 455-459, 1994 ; la méthode des éléments finis (FEM, pour Finite Element Method, en anglais), détaillée dans l'article de Y. YAN, P.L. NUNEZ et R.T. HART, intitulé Finite-element model of the human head: scalp potentials due to dipole sources , Med. Biol. Eng. Comput., vol. 29, no. 5, pp 475-481, 1991; et la méthode des différences finies (FDM pour Finite Difference Method, en anglais), également bien décrite dans l'article de L. LEMIEUX, A. McBRIDE et J.W. HAND, intitulé Calculation of electrical potentials on the surface of a realistic head model by finite differences Phys. Med. Biol., vol. 41, no. 7, pp 1079-1091, 1996. Ajoutons que J. C. MOSHER, R. M. LEAHY, et P. S. LEWIS, dans leurs articles intitulés Matrix kernels for MEG and EEG source localization and imaging , ICASSP 95, 1995 IEEE International Conference on Acoustics Speech and Signal Processing, vol. 5, Detroit, Michigan, mai 1995, pp. 2943û2946 et : EEG and MEG: Forward solutions for inverse methods , IEEE Transactions On Biomedical Engineering, vol. 46, no. 3, pp. 245-259, Mars 1999, ont proposé un mode de calcul commun basé sur i) une formulation matricielle spatio- temporelle du problème direct, et ii) une factorisation de cette formulation matricielle, séparant ainsi les paramètres d'intérêt quasi-linéaires dits de nuisance, à savoir les paramètres d'orientation des sources, des paramètres d'intérêt non linéaires, à savoir les paramètres de localisation desdites sources et des paramètres d'intérêt linéaires, à savoir les décours temporels desdites sources. 2.2 Solutions connues résolvant le problème inverse 2.2.1 Art antérieur connu Quelque soit le domaine d'application, une recherche directe de paramètres multidimensionnels (par exemple de localisation et d'orientation des dipôles) nécessite la résolution d'un problème d'optimisation non convexe. Pour ce faire plusieurs approches ont vu le jour depuis une trentaine d'années. Le domaine des radiocommunications, pour sa part, a et continue d'offrir un panel d'algorithmes traitant ce problème. En télécommunications, les paramètres non linéaires sont les angles d'incidence des sources sur l'antenne de réception et les paramètres de nuisance sont les polarisations des sources qu'il n'est possible d'estimer qu'en présence d'une antenne à diversité de polarisation. Citons les méthodes les plus connues permettant d'estimer les angles d'incidence et les polarisations des sources. L'une des plus célèbres de ces méthodes est la méthode MUSIC (MUltiple Slgnal Classification en anglais, ou classification de signaux multiples, en français) exploitant les statistiques d'ordre 2, décrite par R. O. SCHMIDT dans ses publications : A signal subspace approach to multiple emitter location and spectral estimation , thèse de doctorat de l'Ulniversité de Stanford, Novembre 1981 et Multiple emitter location and signal parameter estimation , IEEE Transactions On Antennas Propagation, vol. 34, no. 3, pp. 276-280, Mars 1986.

Cette méthode s'inscrit au sein d'une grande famille dite des méthodes à sous-espaces qui, en outre, exploitent l'orthogonalité entre l'espace vectoriel des sources et celui du bruit au travers des statistiques d'ordre 2. Ceci est rendu possible à l'ordre 2, lors d'enregistrements multicapteurs, lorsque le nombre de sources, noté P, est strictement plus petit que le nombre d'observations, noté N.

Cependant, bien que la méthode MUSIC de Schmidt permette, en présence d'une antenne à diversité de polarisation, d'estimer les polarisations des sources réceptionnées, elle ne tire pas avantage de la possible factorisation du problème direct. Il faut en effet attendre E. FERRARA et al. dans Direction finding with an array of antennas having diverse polarizations , IEEE Transactions On Antennas Propagation, vol. 31, pp. 231-236, Mars 1983, pour offrir à MUSIC la possibilité de le faire, dissociant ainsi l'estimation des angles d'incidence de celle des polarisations des sources, et de ce fait réduisant considérablement le coût de calcul de l'algorithme. Plusieurs années plus tard, apparaissent des versions séquentielles de MUSIC, exploitant le concept de déflation à l'ordre 2. On compte parmi elles la méthode RapMUSIC, décrite dans l'article de J. C. MOSHER et R. M. LEAHY, Source localization using Recursively Applied and Projected (RAP) MUSIC , IEEE Transactions On Signal Processing, vol. 47, no. 2, pp. 332û340, Février 1999. Cette dernière permet notamment de faciliter la recherche de plusieurs optima locaux dans la métrique MUSIC, en particulier lorsque la dimension de l'espace source augmente. Contrairement aux autres approches séquentielles, la méthode RapMUSIC s'appuie sur les travaux de E. FERRARA, offrant ainsi la possibilité d'exploiter la factorisation de la formulation matricielle du problème direct en présence d'une antenne à diversité de polarisation. Alors que les approches précédentes n'exploitent que les statistiques d'ordre 2, B. PORAT et B. FRIEDLANDER décident d'étendre la version originale de MUSIC aux statistiques d'ordre 4, mais uniquement en contexte non polarisé, tel que décrit dans leur article intitulé Direction finding algorithms based on high-order statistics , IEEE Transactions On Signal Processing, vol. 39, no. 9, pp. 2016-2024, Septembre 1991. La méthode a pour intérêt de permettre la résolution d'un problème inverse mal posé, et plus exactement le traitement d'au plus P = N 2 û 1 sources à partir seulement de N observations. Cette méthode ne permet néanmoins pas, en revanche, d'exploiter une éventuelle factorisation de la formulation matricielle du problème direct.

Par ailleurs, il est important de citer les travaux de M. VIBERG and B. OTTERSTEN, en particulier la technique WSF (Weighted Subspace Fitting) décrite dans leur article intitulé Sensor array processing based on subspace fitting , IEEE Transactions On Signal Processing, vol. 39, pp. 1110-1121, Mai 1991. Il est à noter que cette méthode, dérivée d'une approche de type maximum de vraisemblance, donne des performances, en termes de biais et de variance, approchant en asymptotique celles données par la borne de Cramer-Rao. Quant aux paramètres linéaires, autrement dit les décours temporels des sources, ils peuvent être estimés, uniquement dans le cas surdéterminé, après l'estimation des paramètres non linéaires et de nuisance, notamment en reconstruisant la fonction de transfert liant les observations aux décours temporels des sources. En effet, celle-ci permet, entre autres, de construire le filtre FAS (Filtre Adapté Spatial), décrit dans P. CHEVALIER, Optimal separation of independent narrow-band sources: Concept and Performances , Signal Processing, Elsevier, vol. 73, pp. 27-47, 1999, qui, appliqué aux N observations, permet d'estimer et de reconstruire les décours temporels des P sources. En ce qui concerne la recherche dans le domaine biomédical, plusieurs algorithmes de localisation de sources et de reconstruction d'activités électriques intracérébrales ont vu le jour. Le lecteur pourra se référer à l'article de C. M.

MICHEL, M. M. MURRAY, G. LANTZ, S. GONZALEZ, L. SPINELLI, et R. GRAVE DE PERALTA, EEG source imaging , Clinical Neurophysiology, vol. 115, no. 10, pp. 2195-2222, Octobre 2004, pour une revue des méthodes disponibles. En général, ces algorithmes tentent d'expliquer, d'une manière la plus optimale possible, les potentiels de scalp par les sources intracérébrales. Certaines de ces méthodes permettent de traiter le cas sous-déterminé pour l'estimation à la fois des paramètres linéaires, quasi-linéaires et non linéaires, mais nécessitent cependant l'ajout d'hypothèses a priori sur les sources d'intérêt afin d'obtenir une solution unique au problème. 2.2.2 Inconvénients des solutions connues de l'art antérieur La plupart des méthodes d'estimation de directions d'arrivée sont basées sur les statistiques d'ordre 2 des observations, autrement dit sur les cumulants d'ordre 2 des données acquises au moyen des capteurs utilisés. Ceci implique, de manière explicite ou implicite, de considérer que les sources d'intérêt sont gaussiennes. Or, un inconvénient est qu'une telle hypothèse est très forte, puisque dans beaucoup d'applications les signaux sont généralement non gaussiens et contiennent, en outre, une information statistique pertinente, notamment dans leurs cumulants d'ordre supérieur à 2. Aussi, se contenter d'exploiter les statistiques d'ordre 2 peut de ce fait être limitatif, comme l'ont montré P. CHEVALIER, L. ALBERA, A. FERREOL et P. COMON dans leur article On the virtual array concept for higher order array processing , IEEE Transactions On Signal Processing, vol. 53, no. 4, pp. 1254-1271, Avril 2005, notamment en présence de mélanges sous-déterminés de sources ou bien de bruit gaussien de cohérence spatiale inconnue.

Les méthodes à l'ordre 2 possèdent donc l'inconvénient majeur d'être limitées en termes de performances et ne permettent pas, entre autres, de traiter des mélanges sous-déterminés de sources. En outre, pour ce qui concerne l'algorithme de B. PORAT et B. FRIEDLANDER, s'il offre la possibilité de traiter jusqu'à P = N 2 û 1 sources à partir uniquement de N observations, il ne permet pas, en revanche, de réduire, en tirant par exemple partie d'une éventuelle factorisation de la formulation matricielle du problème direct, le coût de calcul induit par l'optimisation multidimensionnelle. De plus, cette méthode souffre des problèmes justifiant l'emploi des approches séquentielles.

En ce qui concerne la méthode WSF, bien que très performante, les approches de type MUSIC offrent un meilleur compromis entre le niveau de résolution de l'estimation et le coût de calcul engendré. Quant aux méthodes issues du domaine biomédical, et plus particulièrement celles permettant de traiter des mélanges sous-déterminés de sources, elles requièrent en général l'ajout d'hypothèses sur les sources d'intérêt que l'on cherche à étudier. Or ces hypothèses sont parfois purement mathématiques et le cas échéant, souvent déconnectées de la physiologie du problème, ce qui constitue un inconvénient majeur dans l'exploitation et l'interprétation des résultats ainsi obtenus. D'autre part, certaines des méthodes du domaine biomédical nécessitent de reconstruire l'activité électrique en tout point du cerveau susceptible d'être solution du problème inverse, ce qui est très coûteux en termes de calculs, et constitue donc une entrave majeure à la résolution efficace et pertinente de ce dernier. 3. Objectifs de l'invention L'invention a pour objectif de pallier les divers inconvénients de l'état de la technique. Plus précisément, le problème technique que l'invention se propose de résoudre est de fournir une famille de méthodes, baptisées 2q-RapMUSIC (q>_2), permettant simplement et efficacement d'augmenter l'ouverture et d'accroître la résolution d'un réseau de capteurs répartis en certains points d'un environnement multidimensionnel prédéterminé, sans augmenter le nombre de capteurs physiques devant être utilisés, et ce afin de traiter et d'étudier avec une plus grande précision les sources d'intérêt internes à l'environnement considéré et dont le nombre est potentiellement supérieur au nombre de capteurs physiques utilisés, tout en restant insensible à la présence d'un bruit gaussien de cohérence spatiale inconnue. Un tel objectif est particulièrement intéressant dans le domaine biomédical, et plus précisément dans celui de l'électrophysiologie, pour ce qui concerne l'évaluation d'un patient candidat à la chirurgie de l'épilepsie. En effet, l'analyse des données issues de l'imagerie et des observations de surface (EEG et MEG) ne permet pas toujours, au vu des techniques actuelles, de localiser avec précision les zones cérébrales à l'origine des crises d'épilepsie d'un patient, et il est parfois nécessaire d'explorer directement les régions du cerveau impliquées à l'aide d'électrodes intracérébrales. Cependant, l'implantation d'électrodes intracérébrales est un acte chirurgical invasif, donc très contraignant pour le patient. De plus, c'est un acte que seul un nombre limité de médecins ou .2 neurochirurgiens savent pratiquer, ce qui induit des délais d'attente souvent très longs pour les patients. Ainsi, un objectif essentiel de l'invention, appliquée au domaine de l'électrophysiologie, consiste à proposer une technique nouvelle et inventive permettant une localisation fiable et précise sans l'implantation d'électrodes intracérébrales, autrement dit en s'appuyant uniquement sur un réseau de capteurs de surface, dont la pose est simple et totalement non invasive. Une telle technique rendrait également les examens préchirurgicaux exécutables par un plus grand nombre de médecins et permettrait ainsi de réduire les délais d'attente pour les patients. Il est donc bien entendu qu'un autre objectif de l'invention consiste à offrir une telle technique permettant, à partir des seules observations acquises au moyen de capteurs placés en certains points de l'environnement considéré, d'estimer les paramètres linéaires, quasi-linéaires et non linéaires propres aux sources d'intérêt internes audit environnement, de façon à pouvoir surveiller et à pouvoir comprendre aussi bien les comportements normaux de celui-ci (dans le domaine biomédical, l'étude du fonctionnement cérébral dans les phases de sommeil d'un sujet, par exemple), que les comportements anormaux de ce dernier (en sismologie, détection de l'épicentre d'un séisme et surveillance de la propagation de ce dernier, ou dans le domaine biomédical, détection et surveillance d'une crise d'épilepsie, par exemple), sans altération de la qualité des résultats obtenus en comparaison avec les méthodes habituelles, par exemple celles utilisant des électrodes intracérébrales en électrophysiologie. Un objectif supplémentaire de l'invention vise à offrir une telle technique 25 en exploitant une éventuelle factorisation de la formulation matricielle du problème direct afin de réduire sensiblement les coûts de calcul. Encore un objectif de l'invention vise à fournir une telle technique ne nécessitant pas d'identifier les paramètres linéaires des sources en tout point de l'environnement multidimensionnel (par exemple le cerveau) susceptibles d'être 30 solution du problème inverse, de façon à réduire sensiblement le coût de calcul.

Un autre objectif de l'invention est de fournir une telle technique qui permette d'exploiter la théorie des réseaux virtuels d'ordre 2q (q 2) afin de pouvoir i) résoudre un problème mal posé en qui concerne l'identification des paramètres quasi-linéaires et non linéaires, ii) traiter un bruit de nature gaussienne et de cohérence spatiale inconnue, et iii) accroître la résolution de l'estimation des paramètres linéaires, quasi-linéaires et non linéaires en contexte surdéterminé. Un objectif supplémentaire est de fournir une telle technique qui permette d'exploiter le concept de déflation étendu à l'ordre 2q (q ? 2) par les auteurs, afin d'accroître la résolution de l'estimation des paramètres.

Un autre objectif est de fournir une telle invention qui permette de traiter des sources aux paramètres linéaires (décours temporels) potentiellement (mais pas totalement) corrélés. Un dernier objectif de l'invention consiste à fournir une telle technique qui soit d'une part peu coûteuse et d'autre part qui soit relativement simple à mettre en oeuvre, y compris dans des dispositifs, par exemple du type dispositif d'électroencéphalographie, de magnétoencéphalographie, d'étude sismique, ou bien encore de tout type de dispositif dédié à l'étude de l'activité d'un environnement multidimensionnel donné, à partir de données d'observation acquises en certains points de ce dernier. 3. Résumé de l'invention Ces objectifs, ainsi que d'autres qui apparaîtront par la suite, sont atteints selon l'invention à l'aide d'un procédé d'identification de paramètres multidimensionnels, de type linéaires, quasi-linéaires et non linéaires, associés à une pluralité de P?l sources d'intérêt présentes au sein d'un environnement multidimensionnel prédéterminé, au moyen d'une pluralité d'observations en nombre N?1 fini, obtenues à partir de capteurs physiques organisés sous la forme d'un réseau de capteurs répartis en des points prédéfinis dudit environnement.Selon l'invention, les capteurs physiques étant organisés sous la forme d'un réseau de capteurs répartis en des points prédéfinis dudit environnement, le procédé d'identification comprend avantageusement au moins les étapes suivantes : enregistrement de mesures physiques permettant de produire au moins un vecteur de N observations engendrées par un mélange de paramètres linéaires représentatifs de P sources d'intérêt, et d'un bruit additif; construction, à partir dudit au moins un vecteur d'observations, d'une matrice statistique d'ordre 2q (q>_2) des observations ; estimation d'au moins un premier paramètre multidimensionnel desdites P sources d'intérêt au moyen de l'estimée d'au moins un deuxième paramètre multidimensionnel. de façon à analyser et à traiter un plus grand nombre de sources d'intérêt internes audit environnement étudié, à partir d'observations acquises en certains points dudit environnement au moyen d'un nombre plus limité de capteurs physiques. L'invention repose donc sur une approche nouvelle et inventive d'analyse et de traitement de données représentatives de sources d'intérêt internes à un environnement multidimensionnel préalablement défini. Elle constitue un outil de traitement de données puissant, offrant à un utilisateur l'opportunité d'augmenter l'ouverture et d'accroître la résolution de son réseau de capteurs à l'aide de capteurs, dits virtuels, obtenus en exploitant les cumulants d'ordre 2q (q>_2). Tout ceci permet i) d'acquérir une bien meilleure estimation des paramètres d'intérêt notamment en présence de sources partiellement corrélées, ii) d'estimer les paramètres non linéaires et de nuisance en contexte sous-déterminé, et iii) d'être insensible à la présence d'un buit gaussien de cohérence spatiale inconnue. Préférentiellement, le procédé selon l'invention met en oeuvre une étape de localisation et de reconstruction de l'activité électrique générée par la pluralité de P>_1 sources d'intérêt représentatives de sources de courant électrique modélisées sous la forme de dipôles de courant électrique, dits sources dipolaires, lorsque ledit environnement multidimensionnel prédéterminé est un environnement conducteur, les étapes de localisation et de reconstruction tenant compte de la pluralité des N observations en nombre fini. De façon avantageuse, les paramètres linéaires, quasi-linéaires et non linéaires sont respectivement représentatifs des décours temporels ou moments dipolaires, des paramètres d'orientation et de position de chacune des sources de courant électrique. Avantageusement, pour tout instant k, le vecteur d'observations de longueur N s'écrit sous la forme suivante : x(k)= A(0)s(k)+v(k) où : - s(k) est un vecteur, de taille (Pxl), représentatif des paramètres linéaires correspondant aux décours temporels desdites P sources d'intérêt, non gaussiennes et potentiellement (mais pas totalement) corrélées suivant ledit au moins un premier paramètre multidimensionnel ; - A(0) est une matrice de mélange instantané, de taille (NxP), où 10 0={e,,...,epl est l'ensemble des P vecteurs de paramètres quasi- linéaires et non linéaires des sources d'intérêt et où chacun des P vecteurs colonnes de A(61) se décompose sous la forme : a(9) = G(pp) Op avec pp et q représentant respectivement les paramètres non linéaires d'une part et quasi-linéaires d'autre part 15 associés à la p-ième source d'intérêt, la matrice de mélange définissant alors une fonction de transfert entre les P sources d'intérêt et les N observations, et ; v(k) est le vecteur, de taille (Nxl), du bruit additif, indépendant des sources d'intérêt. 20 Préférentiellement, le procédé selon l'invention comporte en outre au moins une étape i) d'estimation des cumulants C,i, e''x d'ordre 2q à partir des K échantillons x(k), ii) de choix du rangement matriciel adéquat pour lequel la matrice statistique d'ordre 2q estimée, de taille (NgxNq), sera notée C2q x . De façon avantageuse, le procédé selon l'invention comporte également au 25 moins une étape d'estimation du rang de la matrice C2g z , et du nombre P de sources impliquées. De façon également avantageuse, le procédé selon l'invention comprend au moins une étape de décomposition en valeurs propres de la matrice C29 x et une étape de construction d'une fonction de coût, dite pseudo-spectre d'ordre 2q ou 30 pseudo-polyspectre, ainsi qu'une étape de minimisation de ladite fonction de coût pour estimer chacun desdits P vecteurs de paramètres quasi-linéaires et non linéaires associé à chacune desdites P de sources d'intérêt, où P est l'estimée de P. Préférentiellement, la fonction de coût s'écrit sous la

det{G' (p)H n9,yG9 (p)} forme J 4 (p) = H , ou : det{Gq (p) G9 (p)} n4,; = (Ê29 ä)H EZq x est un opérateur matriciel, dit projecteur bruit d'ordre 2q, avec Ê29 la matrice des vecteurs propres orthonormalisés associés aux valeurs propres nulles de ladite matrice é;'9 x ; def G'q (p) = G( pr-1 O G(pr' avec G(p) la matrice de gain fonction du vecteur p de paramètres non linéaires et de taille (NxL), où L est la 10 longueur du vecteur de paramètres quasi-linéaires 0.

De façon avantageuse, le procédé selon l'invention met aussi en oeuvre une étape de minimisation de la fonction de coût, réalisée à partir d'un procédé de déflation à l'ordre 2q (q>l) représentatif d'une récursivité d'estimation de chacun des P vecteurs de paramètres quasi-linéaires et non linéaires associé à chacune 15 des P sources d'intérêt.

Préférentiellement, la p-ième étape (1<_p<_P) de la procédure récursive comporte au moins l'une des étapes suivantes :

- recherche du minimum global de la fonction de coût et dont l'estimée est notée 4(p) ;

20 calcul d'un vecteur 4 " (Acp)) en prenant le vecteur propre correspondant à la valeur propre minimale de la matrice t [G,'- A(p H , (P(p))] G9 w (P) )HH4.YG4 p (P(p)) - extraction d'un vecteur 4(p) représentatif d'une estimée du vecteur de paramètres quasi-linéaires 4(p) , à partir du vecteur 4' P' (4(p)) ; 25 - s'il reste au moins une source dont les paramètres quasi-linéaires et non linéaires n'ont pas encore été identifiés, i) construction d'un vecteur def a (B( p)) = G( /;1(p)) 4(p) , puis ii) calcul d'une matrice 4(n)` de taille (N'x1Vq) tenant compte d'un remplacement dudit vecteur a(9(p)) par 17 ledit vecteur a(B(p)) et iii) réaffectation des variables suivant les fonctions suivantes : • P:=Pù1 ; '• 'P' _ q,lo, lo,(P) 9 ( 9,l o )H ~ C lqo,vx e .ù ~q,l oP~ ~ I oPr ~ r 2 ' De façon préférentielle, l'étape d'extraction d'un vecteur les sous-étapes suivantes : - extraction de M=Naù2 vecteurs b~(n) ` (m) de taille (N2xl) ; transformation des vecteurs en M matrices B;(n)` (m) de taille (NxN) ; 10 - calcul d'un vecteur propre commun aux M matrices de l'ensemble Â(r)' et associé à la valeur propre la plus grande. De façon également préférentielle, si le nombre de sources d'intérêt est inférieur au nombre de capteurs physiques, ledit procédé met en oeuvre une étape de construction d'un filtre FAS (pour filtrage Alterné Séquentiel ) défini par la 15 formule par W = [ ê 2 % A (â) , où A (â) est la matrice de mélange reconstruite à partir de l'estimation â desdits paramètres quasi-linéaires et non linéaires des sources d'intérêt, afin d'estimer les paramètres linéaires de ces dernières. L'invention concerne également un dispositif d'identification de paramètres multidimensionnels de sources d'intérêt en ce qu'il comporte un 20 processeur adapté à mettre en oeuvre des étapes du procédé de paramètres multidimensionnels dits linéaires, quasi-linéaires et non linéaires associés à une pluralité de P?1 sources présentes au sein d'un environnement multidimensionnel prédéterminé au moyen d'une pluralité d'observations en nombre N_1 fini, tel que précité. 25 Un tel dispositif peut être intégré, à titre d'exemple illustratif et non limitatif, dans tout type d'appareil d'acquisition de données représentatives de sources d'intérêts de profondeur, à partir d'informations de surface capturées par un ensemble de capteurs, à savoir des appareil du type sismographe dans le domaine de la sismologie, ou bien du type électroencéphalographe et/ou 30 magnétoencéphalographe dans le domaine neurologique. 4(F) comprend L'invention concerne aussi un produit programme d'ordinateur téléchargeable depuis un réseau de communication et/ou stocké sur un support lisible par ordinateur et/ou exécutable par un microprocesseur, un tel programme comprenant des instructions de code de programme pour la mise en oeuvre des étapes du procédé d'identification de paramètres multidimensionnels, de type linéaires, quasi-linéaires et non linéaires, associés à une pluralité de P>_1 sources d'intérêt présentes au sein d'un environnement multidimensionnel prédéterminé, au moyen d'une pluralité d'observations en nombre 1V1 fini, tel que précité. Enfin, l'invention concerne aussi l'application du procédé précité d'identification de paramètres multidimensionnels, du type linéaires, quasis-linéaires et non linéaires, associés à une pluralité de P>_1 sources d'intérêt présentes au sein d'un environnement multidimensionnel prédéterminé, au moyen d'une pluralité d'observations en nombre M_1 fini, aux domaines appartenant aux groupes comprenant : -l'électroencéphalographie ; la magnétoencéphalographie ; - la géophysique ; - la sismologie. Il est bien entendu que l'invention peut s'appliquer à tout autre domaine nécessitant l'identification de paramètres multidimensionnels propres aux sources d'intérêt d'un environnement étudié, à partir uniquement d'informations acquises en certains points de cet environnement, dès lors que le modèle des observations (ou le problème direct) admet une écriture matricielle dissociant les paramètres non linéaires des paramètres quasi-linéaires des dites sources d'intérêt. 5. Liste des figures D'autres caractéristiques et avantages de l'invention apparaîtront plus clairement à la lecture de la description suivante d'un mode de réalisation préférentiel de l'invention, donné à titre d'exemple illustratif et non limitatif, faite en référence aux dessins annexés parmi lesquels : - la figure 1, déjà discutée relativement à l'état de la technique illustre le principe de résolution du problème inverse en électrophysiologie ; la figure 2 illustre le modèle de sources utilisé pour représenter l'activité neuronale ; la figure 3 donne deux exemples de modèles de volume conducteur pouvant être utilisés en EEG et en MEG dans le cadre de l'invention ; la figure 4 présente le principe de l'approximation de Berg et Scherg pour le calcul du potentiel électrique de surface dans le cadre d'un modèle de tête sphérique selon la figure 3. les figures 5 et 6 sont des organigrammes présentant respecteivement les grandes étapes du procédé selon l'invention.

Les différentes figures sont discutées ci-après au travers de la description d'un mode de réalisation détaillée de l'invention. D'ores et déjà, relativement à la figure 2, si l'on considère une synapse excitatrice au niveau des dendrites 20 d'un neurone pyramidal cortical 21, l'activation synaptique provoque au niveau de la membrane postsynaptique une dépolarisation 22 qui peut être assimilée à une entrée de courant. Cette entrée massive d'ions est contrebalancée par des sorties de courant 24 en aval de ce point, le long de la membrane. Un neurone activé peut, par conséquent, être assimilé à un groupe de charges négatives et un groupe de charges positives séparées par une petite distance, c'est à dire à un dipôle de courant. Les courants extracellulaires 24, et par conséquent les champs de potentiels qui s'établissent entre régions positives et négatives, sont à l'origine des activités EEG recueillies en surface. En outre, la figure 3 illustre dans sa partie gauche un modèle de tête sphérique (300) constitué de trois couches concentriques, représentant le cerveau 30, l'os du crâne 31 et la peau du scalp 32. Cette même figure, dans sa partie droite, extraite de la thèse de S. BAILLET (1998), illustre un modèle de tête à géométrie réaliste (301), constitué de 3 milieux (cerveau 33, crâne 34 et scalp 35) et basé, pour chaque patient, sur la segmentation des images IRM anatomique. 6. Description d'un mode de réalisation préféré de l'invention La description du mode de réalisation de l'invention s'inscrit ici dans le domaine de l'EEG, et plus particulièrement de la localisation/reconstruction d'activités électriques intracérébrales à partir des données mesurées à la surface du scalp. Cette description sert de simple exemple illustratif et non limitatif de l'invention. Il est bien entendu possible d'appliquer tout l'enseignement de ce mode de réalisation décrit à tout autre domaine d'application (par exemple à la MEG ou bien à la géophysique) nécessitant à moindre coût de calcul l'estimation de paramètres linéaires, quasi-linéaires et non linéaires propres à des sources d'intérêt internes à un environnement multidimensionnel, à partir de simples données d'observation acquises en certains points dudit environnement, y compris dans des contextes où les sources d'intérêt sont en nombre potentiellement supérieur au nombre des données mesurées, et sont potentiellement (mais pas totalement) cohérentes. 6.1 Hypothèses et modélisation des signaux On observe un ensemble de K vecteurs N-dimensionnels x(k). Chaque vecteur colonne x(k) contient les N différences de potentiels électriques obtenus à l'instant k à partir des N+1 capteurs de surface disposés sur le scalp du patient. On suppose alors que pour tout instant k le vecteur x(k) suit le modèle mathématique suivant : x(k)=A(0)s(k)+v(k) (équation 1) où s(k) est le vecteur des décours temporels des P sources, non gaussiennes et potentiellement (mais pas totalement) corrélées, A(0) est la matrice de mélange instantané de taille (NxP), où 0={8,,...,9) désigne l'ensemble des paramètres quasi-linéaires et non linéaires associés aux P sources, i.e., pour l'application décrite, l'ensemble des paramètres liés à la position et à l'orientation des sources. v(k) est quant-à-lui le vecteur de bruit, supposé gaussien et dont les composantes sont statistiquement indépendantes des décours temporels des sources. Les P sources, comme nous l'avons dit, peuvent être corrélées. On suppose donc qu'elles peuvent être divisées en J groupes, tels que les sources d'un même groupe sont corrélées, alors que les sources de groupes différents sont statistiquement indépendantes. 21 On notera Pj le nombre de sources du j-ème groupe, et P = P.. Notons ~-1 J que le cas où J=P correspond au cas où les P sources sont statistiquement indépendantes alors que le cas où J=1 correspond à celui où elles sont toutes corrélées.

Le vecteur des sources s(k)T s'écrit alors comme la concaténation de J vecteurs sj(k)T. On peut donc écrire que deux composantes de s(k) sont indépendantes si et seulement si elles sont associées à deux vecteurs sj(k) et sj •(k) , tels quel<_j≠j'<_J. La matrice de mélange A(0) s'écrit également comme la concaténation de J matrices AA(Oj), de taille (NxPj), correspondant chacune à un groupe de sources dépendantes. Dans la suite, et par souci de simplicité, les notations O et Oj seront omises. Dans le cadre de la localisation de sources en EEG, chaque vecteur colonne a(8p) de la matrice A(9) peut s'écrire comme le produit d'une matrice de gain G(Op), de taille (Nx3), et d'un vecteur de paramètres quasi-linéaires (ou de nuisance) Op lié à l'orientation de lap-ième source : Vp, 1<_p<_P, a(ee)=G(pp)'p (équation 2) On notera Op = [ppT Op' ]T le vecteur des paramètres non linéaires et quasi-linéaires de la p-ième source, où pp est le vecteur des paramètres non-linéaires de position et 4 le vecteur des paramètres quasi-linéaires d'orientation. Il est important de souligner ici que dès lors qu'une application peut s'appuyer sur le même modèle que celui défini par les équations (1) et (2), les procédés et dispositifs selon l'invention sont applicables. On décrit à présent le modèle de tête sphérique à trois couches, selon la figure 3.a, afin de permettre sa mise en oeuvre. Chaque couche 30, 31, 32 représente un tissu différent, de conductivité Qj (15_ j 3) supposée homogène et isotrope. Toutes les sources de courant sont confinées dans la sphère 30 représentative du cerveau. Elles sont représentées par P dipôles de courant au sens de la figure 2.

Le p-ième dipôle est caractérisé par sa position pp, son orientation Op et son décours temporel sp(k) à l'instant k.

La contribution du p-ième dipôle aux N différences de potentiels électriques de surface s'exprime à travers la matrice de gain G(pp), qui dans le cas d'un modèle de tête sphérique, peut être définie de la façon suivante, en utilisant l'approximation de P. BERG et M. SHERG, présentée dans leur article A fast method for forward computation of multipleshell spherical head models Electroencephalography and Clinical Neurophysiology, vol. 90, no. 1, pp. 58ù64, January 1994, et illustrée dans la figure 4 : 3 3 G(pP)=V ,~ih(ri,f~~PP),K, 2jh(r,N,+,,p;pP) (équation 3) i=l

où V est une matrice dite de switch (selon J. C. MOSHER, R. M. LEAHY, et P. S. LEWIS, Matrix kernels for MEG and EEG source localization and imaging , ICASSP 95, 1995 IEEE International Conference on Acoustics Speech and Signal Processing, vol. 5, Detroit, Michigan, Mai 8-12 1995, pp. 2943ù2946) permettant, dans le cas où les potentiels sont enregistrés par N+l capteurs, de soustraire la valeur recueillie par un des capteurs servant de référence, à l'ensemble des valeurs issues des N autres capteurs, et d'obtenir ainsi N différences de potentiels.

En effet, le potentiel créé par un dipôle dans un milieu à trois sphères homogènes et isotropes peut être approximé par la somme des potentiels créés par trois dipôles, chaque dipôle étant placé dans une seule sphère homogène et isotrope (voir figure 4). Les trois sphères uniques sont identiques à la sphère extérieure (scalp) du modèle à trois couches (même rayon, même conductivité). Les trois dipôles sont orientés de la même façon que le dipôle original , leurs positions et leurs décours temporels sont proportionnels à ceux du dipôle original (coefficients de proportionnalité Âi pour les décours temporels, ,ui pour les positions).

Le vecteur h(r,p), de taille (Nx3) est donné quant-à-lui par l'expression suivante : h(r,p)=[(c,ùc2(r"p))p+czllprri (équation 4) où les paramètres et et c2 sont définis par : 1 (rû p)H p 1 1 = 4~61IPl12,2 + IIrùPIl3 Ilr ù PIl IIrll 2 IIrùPII+IIr11 c2 42r6IIPIl2 ArùPII' + llrllF(r,P) avec a la conductivité de la sphère extérieure (scalp) et : F(r,p) =IIrùPII(IIrlllirùPII+IIriI2 ùPHr) Notons que {À1,.12, /13, P1, p2, p3} constitue l'ensemble des paramètres de Berg définis par P. BERG et M. SHERG, dans leur article intitulé A fast method for forward computation of multiple-shell spherical head models Electroencephalography and Clinical Neurophysiology, vol. 90, no. 1, pp. 58û64, January 1994, et peuvent être obtenus, selon Z. ZHANG, comme décrit dans son article A fast method to compute surface potentials generated by dipoles within multilayer anisotropic spheres , Physics in Medicine and Biology, vol. 40, no. 3, pp. 335û349, Mars 1995, en minimisant la quantité suivante : / \2n-2 ( 3 R' fn û E 2j (1.1 1 û,u; ' ) )2 (équation 7) R3 / j=2 où le nombre de termes N,T,ax doit être suffisamment grand pour garantir une bonne 15 précision de calcul (par exemple N,,,ax = 200), ai et Bi (1 <û 3) sont, respectivement, les conductivités et les rayons des trois sphères, et : Vn, f,= n (équation 8) nm22+(l+n)m21 Les coefficients m21 et m22 sont alors obtenus comme suit : (équation 5) (équation 6) A= Cmu ml2 = 1 2 m21 m22 (2n+ 1)2 k=1n+(n+1)6k/6k+1 (n+1)(6k ' k+1 -1)(R3/Rk)2n+1 n(6k/6k+1ù1)(Rk/R3)2n+1 (n+l)+n(6k/6k+1) où les matrices du produit sont non commutatives. Notons que Mosher et al. ont également donné dans : EEG and MEG: Forward solutions for inverse methods , IEEE Transactions On Biomedical Engineering, vol. 46, no. 3, pp. 245-259, Mars 1999, l'expression de la matrice de gain G(pp) associée au calcul des champs magnétiques dans un modèle de tête sphérique, ainsi qu'une formulation matricielle analogue du problème direct en EEG et en MEG dans le cas des modèles de têtes réalistes (plus exactement pour la méthode BEM ù Boundary Element Method), telles que représentées en figure 3.b. 6.2 Statistiques d'ordre 2q On considère les statistiques d'ordre 2q des observations. On définit, pour tout k, le cumulant d'ordre 2q d'un vecteur x(k) de la façon suivante : c ' " (k) = Cum{x; (k),K ,xq (k),x9 (k),K ,x;Zg (k)'} (équation 10) Cette grandeur statistique est calculée à partir de 2q composantes du vecteur x(k), en considérant q termes xi(k) conjugués et q termes non conjugués. On notera que si les sources d'intérêt et le bruit sont stationnaires, les statistiques d'ordre 2q ne dépendent pas de l'indice de temps k, et pourront être notées C g ,Q x9 . En revanche, si les sources sont cyclostationnaires, cycloergodiques et potentiellement non centrées, d'autres quantités doivent être utilisées telles que celles définies dans l'article de L. ALBERA et al., intitulé Blind Identification of Overcomplete Mixtures of sources (BIOME) , Linear Algebra Applications, Special Issue on Linear Algebra in Signal and Image Processing, vol. 391C, pp. 3-30, November 2004.

Dans la suite, par souci de simplicité, on se limitera à l'étude du cas stationnaire. Le cas non stationnaire peut bien entendu être envisagé au moyen de l'invention.

Il est possible, au moyen du procédé selon l'invention, de ranger l'ensemble des statistiques d'ordre 2q du vecteur x(k) dans une matrice Hermitienne C20 x , de taille (NgxN'), appelée matrice statistique d'ordre 2q. Plusieurs façons de ranger ces statistiques d'ordre 2q dans la matrice C2q x existent. Cependant, seul un nombre fini, noté qo, d'arrangements permet d'obtenir des résultats différents en termes de résolution et de nombre de sources traitées. Ce nombre qo est en général fonction du nombre q précité. Notons cependant que pour des observations à valeurs réelles, telles qu'en EEG, on a qo=l.

Ces qo arrangements sont indexés suivant l'entier -t (051'_<g0û1) et correspondent respectivement à qo matrices statistiques d'ordre 2q, notées CZq.x Le (I, , IZ) -ième élément (1 <_ Nq) de la matrice Czq,x est alors donné par l'expression suivante : Cz'q x (I, ,1'2 1, ) = C'1,K.Îq X (équation 11) où pour chaque nombre tel que OSC ù1 et pour chaque indice if tel que 1<_ü,iz,K'12q N, h =1+EN' i(ij+2q-' -1)+ Nq i-1) i=l (équation 12) Iz =1+EN' i(iJ+q-, -1) + Nq (1 -1) J=l J=1 La propriété de multilinéarité des cumulants et des moments permet d'écrire chaque matrice Czq x (q 1) sous une forme spéciale.

En effet, sous l'hypothèse d'indépendance des sources et du bruit, la matrice `2,.x peut s'écrire, pour chaque 0 51 < (h, sous la forme suivante : CZ'gx =[A q' OA' ']C2 [A q-' Ox A' ' ]H+5(q-l)Cz. (équation 13) où Czq x de taille (NgxNq), C-'q,s de taille (PgxPq) et Czq de taille (NxN) sont respectivement les matrices statistiques d'ordre 2q des vecteurs x(k), s(k) et v(k). 8(.) est le symbole de Kronecker et A q désigne la matrice A élevée à la puissance de Kronecker où O représente l'opérateur du produit de Kronecker. Le nombre e est le même que celui défini dans les équations (11) et (12). On note que l'équation (13) peut être réécrite de la façon suivante : Cza x = E[Aq-' O Aj ~ ~CZe S; [A q-i O A. ' ]H + d(q _ 1)C20a, (équation 14), où la matrice C2q~s de taille (P9 x Pq ), est la matrice statistique d'ordre 2q du vecteur si(k) Il faut noter qu'en pratique les cumulants ne peuvent pas être calculés de manière exacte : il doivent être estimés à partir des composantes du vecteur x(k). En effet, on estime, dans un premier temps, les moments d'ordre m (m<ù2q) à partir des composantes des vecteurs x(k). Puis on utilise la formule de Leonov-Shiryaev décrite pour des données à valeurs réelles dans l'article de P. McCULLAGH, Tensor Methods in Statistics , Chapman and Hall, Monographs on Statistics and Applied Probability, 1987. Cette formule établit la relation mathématique entre les moments d'ordre m (m<_2q) et les cumulants d'ordre 2q. Une formule analogue pour des données à valeurs complexes et pour q appartenant à {2,3} est donnée dans l'article de L. ALBERA et al., intitulé Blind Identification of Overcomplete Mixtures of sources (BIOME) , Linear Algebra Applications, Special Issue on Linear Algebra in Signal and Image Processing, vol. 391C, pp. 3ù30, November 2004. Si les sources d'intérêt et le bruit sont stationnaires et ergodiques, alors les moments d'ordre m (m<ù2q) peuvent être estimés à partir des moyennes temporelles des observations, permettant ainsi d'estimer, via la formule de Leonov-Shiryaev, les cumulants d'ordre 2q. 6.3 Les grandes étapes algorithmiques du procédé selon l'invention Comme l'algorithme 2q-RapMUSIC (q>_2) du procédé selon l'invention, utilisant les statistiques à l'ordre 2q, doit être capable de prendre en compte des signaux potentiellement (mais pas totalement) cohérents, i.e. spatialement corrélés, il est nécessaire de faire les hypothèses suivantes : Al)V15. j<ùJ, P<N; A2) La matrice A q-' O A; ' est de rang plein P9 ; def 1 A3) R(J,q,l) = Ei='Pq <Nq def A4) La matrice A (J,q,1) _ [A, q-' O A'` ` ,K A q ' O A; '] est de rang plein R(J,q,4 En particulier, pour P sources statistiquement dépendantes, c'est-à-dire pour lesquelles J=1, les hypothèses (Al)-(A2) se réduisent donc à : A'1) P<N ; A'2) la matrice A (1, q, l) = A q-' O A' ' est de rang plein Pg.

Alors que pour P sources statistiquement indépendantes, c'est-à-dire pour lesquelles J=P, les hypothèses (Al)-(A2) se réduisent alors à: A"1) P<Nq ; A"2) La matrice A (P, q, l) = A q-' 0A'0' est de rang plein égal à P ; où 0 et A0 représentent respectivement l'opérateur du produit de Khatri-Rao tel que défini dans l'article de R. A. HORN et C. R. JOHNSON intiltulé Topics in Matrix Analysis, Cambridge University Press, New York, 1999 , et A élevé à la puissance q de Khatri-Rao. Dès lors, et tel que décrit ci-après, il devient possible de définir dans le cadre du procédé selon l'invention, le concept de pseudo-spectre d'ordre 2q 15 (q 2 ), ainsi que l'approche dénommée ici par les inventeurs 2q-RapMUSIC . 6.3.1 - Pseudo-spectre d'ordre 2q (q 21 A partir de l'équation (13), il est désormais possible de calculer la décomposition en valeurs propres de la matnce hermitienne C29.x Czq,x =~Ezgs El,e,v][L0 s CEzes Ezgv1H (équation 15) 20 où 4,,s est la matrice réelle de taille (R(J,q,1)xR(J,q,-e)) des valeurs propres non nulles de C' x et EZq s la matrice de taille (N1 xR(J,q,e)) des vecteurs propres orthonormalisés associés. D'autre part, EZq v est la matrice de taille (1Vq x( N'R(J,q,l)) des vecteurs propres orthonormalisés associés aux valeurs propres nulles de Czq,x . Comme la 25 matrice C2q,x est hermitienne, chaque vecteur colonne de EL, est orthogonal à chaque vecteur colonne de EZq,,, . De plus, vect{A (J,q,l )} = vect{EZgs} , par conséquent, chaque vecteur colonne de A (J, q, l) est orthogonal à chaque vecteur colonne de EZq,, .

Ainsi, en notant Bpi le vecteur des paramètres de localisation (position et orientation) de la p-ième (1<ûp<_P) source, appartenant qui plus est au j-ième groupe, et a(B9) q ' Oa(BP)* I le vecteur situé à la "(1ùPjq)(pûl)/(1ûP ) "-ième colonne de la matrice A q-' O A~ ` , les vecteurs a (Bp ) q -1 O a (O" )* I (l et 1 <û j J) sont tous orthogonaux aux vecteurs colonne de Ezq,,, . Ce résultat nous donne lapossibilité de construire une fonction de coût dont la minimisation permet d'estimer les P vecteurs Op' . Cette fonction de coût, baptisée pseudospectre d'ordre 2q (ou pseudo polyspectre), est définie par : J,(0)=[a(B) q-I Oa(9r']Hn ,[a(9) q 1 Oa(9r] (équation 16) où IIq ,, = (Ezq,v )H (EZq ä) est un opérateur matriciel dénommé ici projecteur bruit d'ordre 2q. Il est alors préférable de considérer le critère normalisé suivant J2 (B) = J, (B0l a (B) q-I a(9)' 2 au lieu de celui défini par l'équation 16, afin de rendre le pseudo-polyspectre constant par rapport à o en l'absence de sources. 15 Puis, d'après l'équation (2) et les propriétés du produit de Kronecker, le critère J2 peut être réécrit de la manière suivante : O'HG H G t t l J (B) _ q (p() (p)O 2 0g1HGq (p)Ho' Gq (p)`PgI def iç def q_, ou ~qI =0 q-I 00' et G'(p)=G(p) 0G(p) Minimiser le critère J2 de l'équation (17) par rapport à o revient à i) 20 minimiser la valeur propre minimale, notée aq (p) , de la matrice Gq (p)H I äGq (p) dans la métrique G q (p)" GIq (p) par rapport à p tel que pop, = arg min{2 (p)} , et ii) déduire de Apt le vecteur rhq solution de la p minimisation du critère J2 par rapport à o et notéq' ( pop,) , en prenant le vecteur propre associé à la valeur propre / q(popt) . 25 En effet,

i(popt) et (popt) vérifient : Gq(pop, )H~9,vG4(pop, )P (pop,) ~q( pop, )Gq(pop, )HGq(popt)~gt(pop, ) (équation 18) Par conséquent, la minimisation du critère (17) peut en partie être obtenue en minimisant le critère suivant par rapport à p : I J3 (p) = vpm{[Gq (p)H Gq (p)] Gq (p)H II4,YG9 (p) (équation 19) (équation 17) où vpm{B} désigne la valeur propre minimale de la matrice B. Nous passons ainsi d'une optimisation dans un sous-espace vectoriel dei 6 , à une optimisation dans un sous-espace vectoriel dei 3 , ce qui a pour grand intérêt de réduire le coût de calcul de l'optimisation. Une manière supplémentaire de diminuer le coût de calcul de l'algorithme est de préférer une minimisation du critère suivant à celle de J3 : J4(P) det{G'q(P)HH'q,,G',(P)} det{G9(p)H G9 (p)} où det{B} désigne le déterminant de la matrice B. Remarquons ainsi que de part ces premières étapes du procédé selon l'invention, une recherche multidimensionnelle à la fois en position et en orientation a pu être évitée grâce à l'utilisation du critère .13 défini par l'équation (19). En effet, il suffit désormais d'effectuer une optimisation uniquement par rapport au vecteur de position p sans se soucier de la recherche du vecteur 09 correspondant, puisque ce dernier se déduit ensuite très simplement du vecteur de position optimal trouvé. Puis, une seconde amélioration nous a permis de réduire à nouveau le coût de calcul de l'algorithme en remplaçant la minimisation du critère J3 par celle du critère J4. Effectivement, il est moins coûteux de calculer le déterminant d'une matrice que de décomposer cette dernière en éléments propres. Nous verrons dans la section suivante comment extraire le vecteur d'orientation du vecteur ^P9' pour toutes valeurs de q et de {. 6.3.2 û La méthode 2q-RapMUSIC (q ? 22 On décrit à présent la méthode dénommée 2q-RapMUSIC (q 2) par les inventeurs basée sur i) la factorisation de la formulation matricielle du problème direct, ii) l'exploitation du concept de réseau virtuel d'ordre 2q (q ? 2) au travers du critère J4 défini par l'équation (20), et iii) l'utilisation du concept de déflation étendu à l'ordre 2q et prenant en compte la présence de sources potentiellement corrélées. (équation 20) L'utilisation des ordres supérieurs est nécessaire, comme nous l'avons dit, dans le traitement des mélanges sous-déterminés de sources ou en présence de bruit gaussien de cohérence spatiale inconnue, ou encore afin d'améliorer la résolution de l'estimation des paramètres notamment lorsque les sources sont très proches...DTD: Quant à la nécessité d'utiliser une approche de type déflation, elle vient du fait que le projecteur bruit n'est jamais parfaitement estimé, et que les erreurs d'estimation conduisent nécessairement à rechercher pour les critères JI J4 non pas P minima globaux correspondant respectivement aux P sources, mais plutôt Pù1 minima locaux et un minimum global correspondant par exemple à la source de plus fort Rapport Signal à Bruit (RSB).

Trouver le vecteur position p(,) de la source de plus fort RSB (rapport signal à bruit) est chose facile, et peut être fait en cherchant le minimum de l'équation (20) sur une grille suffisamment échantillonnée des positions susceptibles d'être solutions. Notons que c(.) est un automorphisme de { 1 , 2, ... , P}, autrement dit une permutation de {1, 2, . . . , P}. En effet, les sources ne peuvent être localisées que dans le désordre. Cependant, un coup d'oeil jeté à l'équation (1) suffit pour se rendre compte que modifier l'ordre dans lequel les composantes de s(k) et les colonnes de A(0) correspondantes sont rangées ne change pas l'expression de x(k).

Le vecteur 0q' (pp(,)) , est quant-à-lui obtenu comme le vecteur normalisé qui, multiplié à gauche par la matrice Gq (pp(,)) donne le vecteur q-I .01 a(99(,)) Oa(9,(,)) .

Comme décrit dans la section précédente, ce dernier est obtenu en prenant le vecteur propre associé à la valeur propre minimale de la matrice I [G', (pg(I) )H (pop) )I Gg (e(,) )H nv,vGq (p (I) ) Puis, le vecteur 40) peut être extrait de 4' (pp(,)) . Il faut pour cela extraire en premier lieu les M = N9-2 vecteurs b;(';) (m) (1 m≤ M) de taille (N2x 1) (tels que (p(,)) _ [b;(i) (1)T b;';) (2)TK bop) (M)T ~~ ), puis les convertir en M matrices B) (m) de taille (NxN) (la n-ième colonne de Bq' m) est faite des N éléments consécutifs de b'; )(m) à partir du [N(nû1)+1]-ième), et enfin diagonaliser conjointement l'ensemble 09à,) de matnces définies par {le) (m)BBi) (m)H ,B:i;) (m)T le) (m) , tel que 1<-m<-M, si i'=0, par {B;i;) (m) , tel que 1<_m<-M, si -1)=1 et sinon par {Bj<<) (m)H Bpi) (m),le) (m). B O (m)T , tel que 1<m<-M. Plus précisément, le vecteur 4(,) est estimé par le vecteur propre commun à toutes les matrices de Lgài), et associé à la plus forte valeur propre, est, à un facteur d'échelle près, égal à 40) . Nous parlons d'estimation car le vecteur 4o) ne peut être reconstruit qu'à un facteur d'échelle près, une indétermination intrasèque au problème. Notons cependant que pour la plupart des applications cette indétermination ne constitue pas un problème en soi. En effet, concernant notre problématique de localisation de dipôle de courant, le vecteur 4o) est un vecteur d'orientation. De ce fait, la connaissance du vecteur directeur, de norme unité, de 40) est amplement suffisante, ce qui nous est donné par l'estimation précédente. D'un point de vue algorithmique, la diagonalisation conjointe peut être effectuée dans notre cas de figure au moyen de la méthode JAD (pour Joint Approximate Diagonalization en anglais) proposée par J.-F. CARDOSO et A. SOULOUMIAC dans leur article intitulé Jacobi angles for simultaneous diagonalization , SIAM Journal Matrix Analysis and Applications, vol. 17, no. 1, pp. 161û164, 1996, ou bien de la méthode décrite par A. YEREDOR dans son article intitulé Non-orthogonal joint diagonalization in the least-squares sense with application in blind source separation , IEEE Transactions On Signal Processing, vol. 50, no. 7, pp. 1545û1553, July 2002.

Il est donc facile d'obtenir le vecteur de positions P(,), de même que le vecteur d'orientation correspondant. Il est cependant plus difficile de trouver avec précision les Pû1 minima locaux de J4 restants puisque les techniques usuelles de recherche de minima peuvent manquer des minima superficiels ou ne pas réussir à distinguer deux minima trop voisins.

Ce dernier problème est résolu par le procédé d'identification de paramètres multidimensionnels selon l'invention, par la mise en oeuvre à l'ordre 2q quelconque du concept de déflation, permettant avantageusement, de surcroît, le traitement de sources potentiellement cohérentes.

Concrètement, il s'agit de retirer la contribution de la première source, dont la localisation a été estimée, des observations et de lancer alors une nouvelle recherche de minimum global à partir du pseudo-polyspectre reconstruit à partir des nouvelles observations.

Cette procédure n'en est pas moins très coûteuse car elle nécessite 10 d'estimer P fois la matrice statistique des observations Czq,x .

Une manière moins coûteuse, en termes de temps et de ressources de calcul, consiste à retirer la contribution de la première source estimée non pas des observations mais plutôt de la matrice statistique C2x elle-même.

Cette approche nécessite d'étudier précisément la structure algébrique de 15 la matrice CZq,x , de façon à annuler les colonnes de la matrice A 9-' O A'' clef impliquant le vecteur a(9(1))= G (p(I)) 4(,) à partir de l'équation (13).

Pour se faire, considérons une des propriétés du produit de Kronecker : soient trois vecteurs b, c, d de tailles respectives (Nbxl), (Ncx 1), (Ndx l), le produit de Kronecker vérifie alors bOcOd=(IN, Oc INd)(bOd). 20 En effet, cette propriété permet de montrer la proposition suivante pour deux matrices T'q',70) et Ey,70), de taille (NgxN') définies par, pour tout entier me {0,1,K ,q -1) :

'-W(I) IN' ~I,m-ICI a(9 ), QI q , ) ifmqù1 I TI,m H Ti m m H `~q,el) N9 9,e(1) [(Ti,. 9, I)) 9, (l) (Ti 9,el) ) q,7(I) y rE1q ,me-1) '(I Oa(e(I)) I ) ifm<qù1 N'" , N9-m Alors, multiplier la matrice A q-' O A' ' de taille (NlxPg) à gauche par def E~(1) = ya.f(1)"v. O K ;1'1)) permet d'annuler tous les vecteurs colonnes impliquant a (e,(,)) Le vecteur de position, pf(2) , de la (2)-ième source est alors obtenu en minimisant le critère J4 de l'équation (20), dans laquelle G, p) est désormais remplacé par Eq,G (p) et où II' ., ne sera plus obtenu à partir de la diagonalisation de C2q.X , mais plutôt à partir de la diagonalisation de H 1)C2q,X (Eq,(,) ) H On remarque alors que le rang de la matrice Ee(I)C2qX (E (1) ) est à présent strictement plus petit que R(J,q,f). En effet, nous avons diminué le rang de CZ'q,X en retirant à cette matrice la contribution de la première source estimée. En conséquence, le nombre de vecteurs propres constituant la matrice Ezq,, est augmenté, ce qui traduit une augmentation de la dimension de l'espace bruit, espace engendré par les vecteurs colonnes de E2q.

En rappelant que la matrice Ezq,, intervient directement dans le calcul du projecteur bruit 1110, _ (Ezq )H (E2q,,) , et donc dans le calcul du critère J4 de l'équation (20), nous offrons ainsi la possibilité d'estimer avec une meilleure précision des paramètres de la (2)-ième source. Une fois achevée la minimisation du nouveau pseudo-polyspectre, le 20 vecteur d'orientation de la (2)-ième source, 4(2) , est extrait en suivant la procédure détaillée précédemment pour 4(1) et le vecteur directeur de la c(2)-ième def source, a (ee2 = G(pe2 4(2) peut alors facilement être reconstruit. Il est ensuite procédé par récurrence jusqu'à l'estimation des P vecteurs de paramètres Op = [pET 4 ]T correspondant aux P sources d'intérêt. 25 On précise également que la p-ième étape de la récurrence du procédé selon l'invention consiste à minimiser le critère J4 de l'équation (20) en remplaçant Gq (p) par E (p)K E~~Z~E ~;~Gq (p) et où 11g,, est construit à partir des vecteurs propres associés aux valeurs propres nulles de la 1 q ql H matrice ~~(P)K El 2)Eql g1)C2ql,x(EVp)K El 2)Eg;()) 6.3.3 ù Identifiabilité On peut montrer que le problème du traitement à l'ordre 2q de P sources non gaussiennes potentiellement (mais pas totalement) cohérentes à partir d'un réseau de capteurs générant N observations, est, pour le rangement indexé par e, C21q,x , similaire à un problème de traitement au second ordre desdites P sources à partir d'un réseau virtuel de N q capteurs virtuels dont en général N2eq sont différents. Il est important de noter que le nombre N dépend implicitement de la structure algébrique d u vecteur a(6')= G (p) O , et plus exactement de celle de la matrice de gain G(p) propre à chaque application. De ce fait, pour des valeurs de q et de fixées, le nombre N:q changera d'une application à l'autre. On déduit des résultats précédents que le nombre maximal de sources non gaussiennes et statistiquement indépendantes, i.e. J=P, pouvant être traitées par l'algorithme 2qûRapMUSIC pour le rangement indexé par est9 1. 6.4 Résumé des grandes étapes successives du procédé selon l'invention Relativement à la figure 5, on résume cidessous les différentes étapes du procédé selon l'invention en supposant que K échantillons temporels d'observations de surface, x(k) (1 k K ), sont disponibles. Etape 1 (51) : Choisir l'ordre statistique 2q approprié en fonction du nombre P de sources que l'on souhaite traiter. En pratique, q est la valeur minimale qui assure le traitement de toutes les sources potentiellement présentes dans l'environnement multidimensionnel étudié.

Etape 2 (52) : Estimer les cumulants C; `''x d'ordre 2q à partir des K échantillons x(k) et choisir le rangement matriciel adéquat pour lequel la matrice statistique d'ordre 2q estimée sera notée C21q x Etape 3 (53) : Calculer les valeurs propres de la matrice hermitienne Czq X et extraire une estimée Ê29 , de la matrice Ezq . Cette étape peut nécessiter l'estimation du rang de la matrice C''2;''; pour des cas de figure où le nombre de sources et/ou leur cohérence spatiale sont/est inconnu(e)(s) (dès lors nous noterons P l'estimée de P).

Etape 4 (54) : Calculer une estimée, 119,v = (ÊZq x )H E29 x , du projecteur bruit H , d'ordre 2q.

Etape 5 (55) : Calculer une estimée, du critère J4 de l'équation 20 en utilisant la matrice fi , pour une grille de positions appropriée selon les termes

précités, et chercher son minimum global, noté .

Etape 6 (56) : Calculer le vecteur 4' ` (4(,)) en prenant le vecteur propre correspondant a la valeur propre minimale de la matrice 1 [G,''P 40 H G q (4(i) )] G' q (4(,) )H HG' q' (Pf(i) )

Etape 7 (57) : Extraire le vecteur 4(,), estimée du vecteur quasi-linéaire 4(1), à partir du vecteur îq' (4(1)) . Pour cela, extraire tout d'abord les M = Nq-2 vecteurs b;(';) ` (m) de taille (N2x1), puis les transformer alors en M matrices Ê"" (m) de taille (NxN), et enfin calculer le vecteur propre commun aux M matrices de l'ensemble et associé à la valeur propre la plus grande. Utiliser pour cela un algorithme de diagonalisation conjointe de matrices tels que ceux précités.

Etape 8.1 (581) : Si les P vecteurs de paramètres quasi-linéaires et non linéaires des sources ne sont pas encore tous identifiés, construire le vecteur def a(O (1)) = G(4(1))4(,) et calculer alors la matrice de taille (Ngx1Vq) de la manière décrite dans la section précédente en remplaçant le vecteur a (Be(1)) par a (B.(,) ) Etape 8.2 (582) : Puis, réaffecter les variables de la manière suivante : P:=Pû1 ; - GQ P`(P):=Ê(q`G'9 (P)

- ê1 ,' tq,1-,c1 P Êq,1- H

2q,x ' ~(1) 2q,x ( f(1) et retourner à l'étape 3 (53), Sinon aller à l'étape 9 (59).25 Etape 9 (59): Si P≤. N (c'est-à-dire si le mélange de sources est surdéterminé), estimer les P paramètres linéaires sp(k) pour chaque valeur k en appliquant au vecteur d'observations x(k) le filtre FAS défini par à(' =~â 1-1 A(ê), où A(ê) est la matrice de mélange reconstruite à partir de 5 l'estimation ê des paramètres quasi-linéaires et non linéaires des sources en utilisant l'équation (2).

Claims (13)

REVENDICATIONS

1. Procédé d'identification de paramètres multidimensionnels, de type linéaires, quasi-linéaires et non linéaires, associés à une pluralité de P>_1 sources d'intérêt présentes au sein d'un environnement multidimensionnel prédéterminé, au moyen d'une pluralité d'observations en nombre N?l fini, obtenues à partir de capteurs physiques organisés sous la forme d'un réseau de capteurs répartis en des points prédéfinis dudit environnement, caractérisé en ce que ledit procédé d'identification comprend au moins les étapes suivantes : enregistrement de mesures physiques permettant de produire au moins un vecteur de N observations engendrées par un mélange de paramètres linéaires représentatifs de P sources d'intérêt, et d'un bruit additif; - construction, à partir dudit au moins un vecteur d'observations, d'une matrice statistique d'ordre 2q (q>_2) des observations ; - estimation d'au moins un premier paramètre multidimensionnel desdites P sources d'intérêt, au moyen de l'estimée d'au moins un deuxième paramètre multidimensionnel.

2. Procédé d'identification de paramètres multidimensionnels de sources d'intérêt selon la revendication 1, caractérisé en ce qu'il met en oeuvre une étape de localisation et de reconstruction de l'activité électrique générée par ladite pluralité de P?1 sources d'intérêt représentatives de sources de courant électrique modélisées sous la forme de dipôles de courant électrique, dits sources dipolaires, lorsque ledit environnement multidimensionnel prédéterminé est un volume conducteur, et en ce que lesdites étapes de localisation et de reconstruction tiennent compte de la pluralité desdites N observations en nombre fini.

3. Procédé d'identification de paramètres multidimensionnels de sources d'intérêt selon la revendication 2, caractérisé en ce que les paramètreslinéaires, quasi-linéaires et non linéaires sont respectivement représentatifs des décours temporels ou moments dipolaires, des paramètres d'orientation et de position de chacune desdites sources de courant électrique.

4. Procédé d'identification de paramètres multidimensionnels de sources d'intérêt selon l'une quelconque des revendications 1 à 3, caractérisé en ce, pour tout instant k, le vecteur d'observations de longueur N s'écrit sous la forme suivante : x(k)=A(0)s(k)+v(k) où : - s(k) est un vecteur, de taille (Pxl), représentatif des paramètres linéaires correspondant aux décours temporels desdites P sources d'intérêt, non gaussiennes et potentiellement (mais pas totalement) corrélées suivant ledit au moins un premier paramètre multidimensionnel ; - A(C est une matrice de mélange instantané, de taille (NxP), 15 où 0=10,,...41 est l'ensemble des P vecteurs de paramètres quasi-linéaires et non linéaires desdites sources d'intérêt, et où chacun des P vecteurs colonnes de A(6 se décompose sous la forme : a (Op) = G( pp ) ~p , avec pp et 4, représentant respectivement les paramètres non linéaires et 20 quasi-linéaires associés à la p-ième source d'intérêt, ladite matrice de mélange définissant une fonction de transfert entre lesdites P sources d'intérêt et lesdites N observations, et ; - v(k) est le vecteur, de taille (Nx 1), du bruit additif, 25 indépendant desdites sources d'intérêt.

5. Procédé d'identification de paramètres multidimensionnels de sources d'intérêt selon l'une quelconque des revendications 1 à 4, caractérisé en ce qu'il comporte au moins une étape i) d'estimation des cumulants d'ordre 2q à partir des K échantillons x(k), ii) de choix du 30 rangement matriciel adéquat pour lequel la matrice statistique d'ordre 2q estimée, de taille (Ni' xN4), sera notée Czgx .

6. Procédé d'identification de paramètres multidimensionnels de sources d'intérêt selon l'une quelconque des revendications 1 à 5, caractérisé en ce qu'il comporte au moins une étape d'estimation du rang de ladite estimée Czg X ,et du nombre P de sources impliquées.

7. Procédé d'identification de paramètres multidimensionnels de sources d'intérêt selon la revendication 6, caractérisé en ce qu'il comprend au moins une étape de décomposition en valeurs propres de ladite estimée C29 x et une étape de construction d'une fonction de coût, dite pseudo-spectre d'ordre 2q ou pseudo-polyspectre, et une étape de minimisation de ladite fonction de coût pour estimer chacun desdits P vecteurs de paramètres quasi-linéaires et non linéaires associés à chacune desdites P de sources d'intérêt, où P où est l'estimée de P.

8. Procédé d'identification de paramètres multidimensionnels de sources d'intérêt selon la revendication 7, caractérisé en ce que ladite fonction det{G` ( p ) H II' G` ( p) } de coût s'écrit sous la forme 1'4( p) = q'v q où : det{G`q ( p ) H Gq (p)} - ui ; = ,-,,x )H .4 X est un opérateur matriciel, dit projecteur bruit d'ordre 2q, avec ÊY la matrice des vecteurs propres orthonormalisés associés aux valeurs propres nulles de ladite matrice C2l9 X ; - G(p)=G(p) O G(p)' ' avec G(p) la matrice de gain fonction du vecteur p de paramètres non linéaires et de taille (NxL), où L est la longueur du vecteur de paramètres quasi-linéaires 0.

9. Procédé d'identification de paramètres multidimensionnels de sources d'intérêt selon l'une quelconque des revendications 7 et 8, caractérisé en ce qu'il met en oeuvre une étape de minimisation de ladite fonction de coût, réalisée à partir d'une déflation à l'ordre 2q (q>1) représentative d'une récursivité d'estimation desdits P vecteurs def qù1contenant lesdits paramètres quasi-linéaires et non linéaires associés à chacune desdites P sources d'intérêt.

10. Procédé d'identification de paramètres multidimensionnels de sources d'intérêt selon la revendication 9, caractérisé en ce que la p-ième étape (1<ùp<_P) de la procédure récursive comporte au moins l'une des étapes suivantes: - recherche du minimum global de ladite fonction de coût et dont l'estimée est notée p(P) ; - calcul d'un vecteur ' r (peP en prenant le vecteur propre correspondant à la valeur propre minimale de la matrice [Gq'm' fI(p H G4.,r (P~(P) )~ G q r (p (P) )H Hg YGep (2 (P) ) -extraction d'un vecteur 4(P) représentatif d'une estimée du vecteur de paramètres de nuisance 4(P) , à partir du vecteur 49 .Pr 4(P) s'il reste au moins une source dont les paramètres quasi-linéaires et non linéaires n'ont pas encore été identifiés, i) def construction d'un vecteur a(B(P)) = G(peP))4(P) , puis ii) calcul d'une matrice 4(n) de taille (1VgxNq) tenant compte d'un remplacement dudit vecteur a(O(P)) par ledit vecteur 20 a (B(P)) , et iii) réaffectation des variables suivant les fonctions suivantes : • P:=Pù1 ; • e ` (P) = sepj `G9 `r (P) C29+x := iavr -4,1orrC1ow ~9+iopr H • (1) 24,x ( 1) ) 25

11. Procédé d'identification de paramètres multidimensionnels de sources d'intérêt selon la revendication 10, caractérisé en ce que ladite étape d'extraction d'un vecteur 4(P) comprend les sous-étapes suivantes : -extraction de M= Nqù2 vecteurs b;(P)r (m) de taille (N2x 1) ; transformation desdits vecteurs en M matrices b';(' pl lm) de 30 taille (NxN) ; 10 15

12. 513. 14. 15. 30- calcul d'un vecteur propre commun aux M matrices de l'ensemble A(P)' et associé à la valeur propre la plus grande. Procédé d'identification de paramètres multidimensionnels de sources d'intérêt selon l'une quelconque des revendications 1 à 11, caractérisé en ce que, lesdites P sources d'intérêt sont inférieures en nombre aux observations, il met en oeuvre une étape de construction d'un filtre FAS (pour filtrage Alterné Séquentiel ) défini par la formule par W [ê% Tl A(â) où A(â) est la matrice de mélange reconstruite à partir de l'estimation â desdits paramètres quasi-linéaires et non linéaires desdites sources d'intérêt et d'estimation des paramètres linéaires de ces dernières. Dispositif d'identification de paramètres multidimensionnels de sources d'intérêt en ce qu'il comporte un processeur adapté à mettre en oeuvre des étapes du procédé de paramètres multidimensionnels dits linéaires, quasi-linéaires et non linéaires associés à une pluralité de P?1 sources présentes au sein d'un environnement multidimensionnel prédéterminé, au moyen d'une pluralité d'observations en nombre N?1 fini, selon l'une au moins des revendications 1 à 12. Appareil d'électroencéphalographie et/ou de magnétoencéphalographie, caractérisé en ce qu'il comprend un dispositif d'identification de paramètres multidimensionnels de sources d'intérêt selon la revendication

13. Produit programme d'ordinateur téléchargeable depuis un réseau de communication et/ou stocké sur un support lisible par ordinateur et/ou exécutable par un microprocesseur, caractérisé en ce qu'il comprend des instructions de code de programme pour la mise en oeuvre des étapes du procédé d'identification de paramètres multidimensionnels, de type linéaires, quasi-linéaires et non linéaires, associés à une pluralité de P?l sources d'intérêt présentes au sein d'un environnement multidimensionnel prédéterminé, au moyen d'une 10 15pluralité d'observations en nombre M_1 fini, selon l'une quelconque des revendications 1 à 12. 16. Application du procédé d'identification de paramètres multidimensionnels, de type linéaires, quasi-linéaires et non linéaires, associés à une pluralité de P?1 sources d'intérêt présentes au sein d'un environnement multidimensionnel prédéterminé, au moyen d'une pluralité d'observations en nombre N?1 fini, selon l'une quelconque des revendications 1 à 12, aux domaines appartenant aux groupes comprenant : -l'électroencéphalographie ; - la magnétoencéphalographie ; la géophysique ; la sismologie. 20