WO2007057453A1 - Procede et dispositif d'identification de parametres multidimensionnels : application a la localisation et la reconstruction d'activites electriques de profondeur au moyen d'observations de surface - Google Patents

Procede et dispositif d'identification de parametres multidimensionnels : application a la localisation et la reconstruction d'activites electriques de profondeur au moyen d'observations de surface Download PDF

Info

Publication number
WO2007057453A1
WO2007057453A1 PCT/EP2006/068636 EP2006068636W WO2007057453A1 WO 2007057453 A1 WO2007057453 A1 WO 2007057453A1 EP 2006068636 W EP2006068636 W EP 2006068636W WO 2007057453 A1 WO2007057453 A1 WO 2007057453A1
Authority
WO
WIPO (PCT)
Prior art keywords
sources
parameters
linear
interest
matrix
Prior art date
Application number
PCT/EP2006/068636
Other languages
English (en)
Inventor
Laurent Albera
Delphine Cosandier-Rimele
Isabelle Merlet
Fabrice Wendling
Original Assignee
Universite De Rennes 1
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Universite De Rennes 1 filed Critical Universite De Rennes 1
Priority to US12/094,078 priority Critical patent/US20090093964A1/en
Publication of WO2007057453A1 publication Critical patent/WO2007057453A1/fr

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/22—Source localisation; Inverse modelling

Definitions

  • Method and device for identifying multidimensional parameters application to the location and reconstruction of electrical depth activities by means of surface observations
  • the field of the invention is that of the acquisition and processing of signals representative of activities generated by a set of sources internal to a given multidimensional environment that one seeks to study.
  • the invention relates to a new multidimensional parameter identification technique that allows, among other things, the location and reconstruction of electrical activities, commonly called sources, generated within a multidimensional environment, from observations alone. acquired at certain points of said environment by means of a set of physical sensors.
  • the invention is therefore part of a context and a technical problem that may interest many disciplines for which this problem of identifying multidimensional parameters of sources of interest from observations is essential and necessary for the better understanding of internal phenomena. to the studied environment.
  • the EEG collects surface electrical activity by means of sensors (recording electrodes) arranged in a standardized manner on the scalp and
  • RECTIFIED SHEET (RULE 91) ISA / EP allows to represent the evolution over time of a difference of electrical potentials between each electrode and a common reference.
  • the MEG records the magnetic field (of very low amplitude) induced by cerebral electrical activity thanks to ultrasensitive sensors placed on a helmet covering the entire scalp.
  • these two techniques make it possible to record particular spontaneous rhythms linked to normal physiological states (sleep, sleep, attention); they also make it possible to collect the responses evoked by a given stimulation, which reflect the activation of the different brain structures involved in perceptual and / or cognitive processes.
  • the EEG and MEG recordings make it possible to better understand the mechanisms involved in the context of certain neurological pathologies such as epilepsy, Alzheimer's disease or brain tumors.
  • the analysis of the neural bases of a normal complex cognitive system or of a pathological dysfunction, such as an epileptic seizure, for example requires the precise identification of both the participating brain regions and the temporal sequence of activation between these regions.
  • the known techniques of functional imaging fMRI for functional magnetic resonance imaging, PET for Positron Emission Tomography, TEMP for Mono Photonic Emission Tomography
  • these methods remain severely limited in their ability to reveal the temporal information of this activation since, at best, they provide an average activation image over several hundred milliseconds.
  • RECTIFIED SHEET (RULE 91) ISA / EP requires the use of mathematical tools that are both reliable and accurate to solve the opposite problem, that is to say to locate and reconstruct the sources of brain activity considered from only surface observations.
  • the inverse problem in human electrophysiology is defined as the possibility, from simple surface recordings (electrical potentials and / or magnetic fields) 11 of brain activity derived from sensors 10 and using models 12 and sources appropriate to the conduction medium considered, to identify the brain regions 13 responsible for the recorded EEG and / or MEG activities. More specifically, the inverse problem in EEG and MEG consists of estimating the parameters of the dipolar sources of brain activity and reconstructing the associated temporal paths, from surface observations alone.
  • the most commonly used source model is the dipolar model, which assimilates the activity of a small cortical zone to that of a current dipole, as shown in FIG. distinguishes then two cases.
  • the number of dipoles is assumed to be less than or equal to the number of surface observations (EEG or MEG activities). This is referred to as an overdetermined mixture of sources, or a "well-posed" problem. However, it is inaccurate in practice to reduce normal or pathological brain function to the activity of a limited number of sources.
  • the electric potential and the magnetic field recorded on the surface of the scalp of a subject are also related to the physical and geometric constraints of the various tissues of the head.
  • ISA / EP is assimilated to a conductive volume that must take into account the inhomogeneities of the different media (brain, cerebrospinal fluid, skull and skin).
  • the head is thus generally modeled by a set of three or four concentric layers, of different conductivities, representing the different tissues crossed when the signal from the sources reaches the skin surface.
  • each of these media can be considered as isotropic, as described in the article by S. RUSH and D. A. DRISCOLL, "EEG electrode sensitivity - An application of reciprocity," IEEE Transactions On Biomedical Engineering, vol. 38, pp. 15-22, January 1969, or else anisotropic, as described in the article by J. C. De MUNCK and M. J. PETERS,
  • the simplest model of geometry is the spherical model (figure 3.a) which assimilates the head to a set of concentric spheres, where each layer corresponds to a different fabric.
  • the most commonly used model is a three-sphere model which represents respectively the brain 30, the skull 31 and the scalp 32 of a subject.
  • the spherical model is, however, only a rough approximation of the geometry of the head. Models with a more realistic geometry (figure 3.b) have therefore been developed, and are constructed, for each subject, from anatomical MRI images. Methods make it possible to extract, by segmentation of the MRI images, the contours of the three structures of interest that are the brain 33, the skull 34 and the scalp 35, and then to generate 3D meshes of these three surfaces. 2.1.3 Solving the direct problem
  • FDM finite difference method
  • the nonlinear parameters are the angles of incidence of the sources on the receiving antenna and the nuisance parameters are the polarizations of the sources which can only be estimated in the presence of a polarization diversity.
  • This method is part of a large family called sub-space methods which, in addition, exploit the orthogonality between the vector space of the sources and that of the noise through the statistics of order 2. This is made possible at order 2, during multisensor recordings, when the number of sources, noted P, is strictly smaller than the number of observations, noted N.
  • the MUSIC method of Schmidt allows in the presence of a polarization diversity antenna, to estimate the polarizations of the received sources it does not take advantage of the possible factoring of the direct problem.
  • the MUSIC method makes it possible to estimate both the nonlinear parameters and the quasi-linear parameters, this is difficult to achieve in the operational context due to the computational complexity that this entails.
  • E. FERRARA proposes a new version of MUSIC that exploits the factorization of the matrix formulation of the direct problem, and dissociates the estimation of the nonlinear parameters (that is to say the positions of the current dipoles in human electrophysiology) of the quasi-linear parameters (that is to say the orientations of said dipoles) of the sources, thereby considerably reducing the cost of calculating the algorithm.
  • the algorithm of E. FERRARA then makes feasible in practical context the simultaneous estimation of nonlinear and quasi-linear parameters.
  • the RapMUSIC method is based on the work of E. FERRARA, thus offering the possibility of exploiting a possible factorization of the matrix formulation of the direct problem in the presence of a polarization diversity antenna, for example , and more generally when the quasi-linear parameters are unknown and must be estimated.
  • this method does not make it possible to exploit a possible factorization of the matrix formulation of the direct problem and is algorithmically different from the original version of MUSIC not only by the use of higher order statistics, but also by by the fact that the algebraic structure of the quadricovariance matrix is different from that of the covariance matrix, and then asks the authors to modify the original version of MUSIC.
  • FERRARA's MUSIC algorithm or even the RapMUSIC method without modifying the latter when a single contracted quadricovariance matrix is used.
  • the authors propose to estimate the signal space and the noise space by joint diagonalization of the different matrices of contracted quadricovariance used.
  • the MUSIC and RapMUSIC algorithms can then be used again without the need to modify them.
  • the linear parameters in other words the temporal decours of the sources, they can be estimated, only in the overdetermined case, after the estimation of the nonlinear, quasi-linear and nuisance parameters, in particular by reconstructing the transfer function linking the observations to the temporal course of sources. Indeed, it allows, among other things, to build the filter FAS (Adapted Spatial Filter), described in P. CHEVALIER,
  • SATOSHI NIUIMA and SHOOGO UENO assume in the implementation of their method that the sources are spatially statistically independent to order 4, which in operational context is far from always being verified.
  • two epileptic zones may have highly correlated electrical activities and yet be several centimeters apart.
  • the MUSIC-type approaches offer a better compromise between the resolution level of the estimate and the calculation cost generated.
  • the object of the invention is to overcome the various disadvantages of the state of the art.
  • the technical problem that the invention proposes to solve is to provide a family of methods, called "2g-RapMUSIC"
  • RECTIFIED SHEET (RULE 91) ISA / EP predetermined multi-dimensional, without increasing the number of physical sensors to be used, and in order to process and study with greater precision the sources of interest internal to the environment and whose number is potentially greater than the number of sensors used, while remaining insensitive to the presence of a Gaussian noise of unknown spatial coherence.
  • Such an objective is particularly interesting in the biomedical field, and more specifically in the field of electrophysiology, as regards the evaluation of a patient candidate for epilepsy surgery.
  • the analysis of data from imaging and surface observations does not always make it possible, according to current techniques, to precisely locate the brain areas at the origin of epileptic seizures. of a patient, and it is sometimes necessary to directly explore the regions of the brain involved using intracerebral electrodes.
  • implantation of intracerebral electrodes is an invasive surgical procedure, therefore very restrictive for the patient.
  • an essential objective of the invention applied to the field of electrophysiology, consists in proposing a new and inventive technique allowing a reliable and precise localization without implantation of intracerebral electrodes, in other words relying solely on a network of surface sensors, whose installation is simple and completely non-invasive.
  • Such a technique would also make the pre-surgical examinations more feasible for a larger number of physicians and thus reduce wait times for patients.
  • Another objective of the invention is to offer such a technique that, from the observations obtained by means of sensors placed at certain points of the environment under consideration, to estimate the linear parameters, quasi -linear and nonlinear specific to the internal sources of interest in said environment, so as to be able to monitor and be able to understand both the normal behaviors of the environment (in the field
  • RECTIFIED SHEET (RULE 91) ISA / EP biomedical, the study of brain function in the sleep phases of a subject, for example), that the abnormal behaviors of the latter (in seismology, detection of the epicenter of an earthquake and monitoring of the propagation of the latter , or in the biomedical field, detection and monitoring of an epileptic seizure, for example), without altering the quality of the results obtained in comparison with the usual methods, for example those using intracerebral electrodes in electrophysiology.
  • Another objective of the invention is to provide such a technique that makes it possible to exploit the theory of virtual networks of order 2q (q ⁇ 2), and / or, unlike the UFO-MUSIC method, the algebraic structure of the tensor cumulative order 2q (q ⁇ 2) in its entirety and not that of matrix slices of said tensor, in order to simultaneously i) solve an ill-posed problem with regard to the identification of quasi-linear and non-linear parameters, ii) treat noise of a Gaussian nature and unknown spatial coherence, and iii) increase the resolution of the estimation of linear, quasi-linear and non-linear parameters in an overdetermined context.
  • a further object of the invention is to provide such a technique capable of exploiting a possible factorization of the matrix formulation of the direct problem in order to dissociate the estimation of the nonlinear parameters from the sources of the estimate of the quasi-linear parameters, and to reduce significantly the calculation costs when the quasi-linear parameters are unknown.
  • This objective aims to guarantee the feasibility of the process in operational context when the quasi-linear parameters of the sources are unknown.
  • An additional objective is to provide such a technique that allows to exploit the concept of extended deflation and formalized to the order 2q (q ⁇ 2) by the authors, in order to increase the resolution of the estimation of the parameters.
  • Another objective is to provide such an invention that makes it possible to process sources with linear parameters (time course) potentially (but not totally) correlated, especially in space.
  • Yet another object of the invention is to provide such a technique not exploiting any other information a priori on the linear parameters of the sources.
  • a last objective of the invention is to provide such a technique which is on the one hand inexpensive and on the other hand which is relatively simple to implement, including in devices, for example of the electroencephalographic device type, magnetoencephalography, seismic study or exploration, or even any type of device dedicated to the study of the activity of a given multidimensional environment, from observation data acquired at certain points of the latter .
  • the identification method advantageously comprises at least the following steps:
  • RECTIFIED FEUJLLE (RULE 91) ISA / EP estimating at least a first multidimensional parameter of said P sources of interest by means of the estimate of at least one second multidimensional parameter. in order to analyze and process a greater number of internal sources of interest in said environment studied, based on observations acquired at certain points of said environment by means of a more limited number of physical sensors.
  • the invention is therefore based on a novel and inventive approach for analyzing and processing data representative of internal sources of interest in a previously defined multidimensional environment. It is a powerful data processing tool, offering a user the opportunity to increase the opening and increase the resolution of its sensor network using so-called virtual sensors, obtained by exploiting the cumulative data. order 2q (q ⁇ 2). All this allows i) to acquire a much better estimate of the parameters of interest especially in the presence of partially correlated sources, ii) to estimate the nonlinear and nuisance parameters in an under-determined context, and iii) to be insensitive to the presence of a Gaussian buit of unknown spatial coherence.
  • the method according to the invention implements a step of localization and reconstruction of the electrical activity generated by the plurality of P ⁇ i sources of interest representative of electric current sources modeled in the form of electric current dipoles.
  • said dipolar sources when said predetermined multidimensional environment is a conductive environment, the location and reconstruction steps taking into account the plurality of N observations in finite number.
  • the linear, quasi-linear and non-linear parameters are respectively representative of the time steps or dipole moments, of the orientation and position parameters of each of the electric current sources.
  • - v (k) is the vector, of size (NxI), of the additive noise, independent of the sources of interest.
  • the method according to the invention also comprises at least one step i) of estimating cumulants C 1 "* 1 ";”; " Of order 2q from K samples x (k), ii) of choice of adequate matrix storage for which the statistical matrix of order 2q estimated, of size (WxN * 1 ), will be denoted by C ⁇ x .
  • the method according to the invention also comprises at least one step of estimating the rank of the matrix C ⁇ x , and the number P of sources involved.
  • the method according to the invention comprises at least one step of decomposition into eigenvalues of the matrix C ⁇ x and a step of construction of a cost function, called pseudo-spectrum of order 2q or pseudo-spectrum of polyspectrum, and a step of minimizing said cost function to estimate each of said P quasi-linear and non-linear parameter vectors associated with each of said P of sources of interest, where P is the estimate of P.
  • the method according to the invention also implements a step of minimizing the cost function, carried out from a 2q (q> 1) deflation process representative of a estimation recursion of each of the P quasi-linear and non-linear parameter vectors associated with each of the P sources of interest.
  • the /? - th step (l ⁇ p ⁇ P) of the recursive procedure comprises at least one of the following steps:
  • the step of extracting a vector 4 ⁇ comprises the following substeps:
  • FAS Simple Alternate Filtering
  • the invention also relates to a device for identifying multidimensional parameters of sources of interest in that it comprises a processor adapted to implement steps of the method of multidimensional parameters called linear, quasi-linear and non-linear associated with a plurality of P ⁇ sources present within a predetermined multidimensional environment by a plurality of N> 1 finite number observations, as above.
  • a device can be integrated, by way of illustrative and nonlimiting example, in any type of data acquisition apparatus representative of sources of depth interest, from surface information captured by a set of sensors. that is, seismograph-type apparatus in the field of seismology, or electroencephalograph and / or magnetoencephalograph in the neurological field.
  • the invention also relates to a computer program product downloadable from a communication network and / or stored on a computer readable medium and / or executable by a microprocessor, such a program comprising program code instructions for the implementation steps of the method for identifying multidimensional parameters of the linear, quasi-linear and non-linear type associated with a plurality of P ⁇ l sources
  • the invention also relates to the application of the aforementioned method for identifying multidimensional parameters of the linear, quasislinear and nonlinear type associated with a plurality of P ⁇ l sources of interest present in an environment.
  • predetermined multidimensional array by means of a plurality of finite numbered observations N ⁇ 1, to the domains belonging to the groups comprising:
  • FIG. 1, already discussed in relation to the state of the art illustrates the principle of solving the inverse problem in electrophysiology
  • Figure 2 illustrates the source model used to represent neuronal activity
  • FIG. 3 gives two examples of conductive volume models that can be used in EEG and MEG in the context of the invention
  • FIG. 4 presents the principle of the Berg and Scherg approximation for the calculation of the surface electrical potential in the context of a spherical head model according to FIG.
  • FIGS. 5 and 6 are flowcharts presenting the major steps of the process according to the invention in a respectful manner.
  • Figure 3 illustrates in its left part a spherical head model (300) consisting of three concentric layers, representing the brain
  • Each column vector x (k) contains the N electric potential differences obtained at time k from the iV + 1 surface sensors arranged on the scalp of the patient.
  • v (k) is the noise vector, supposedly Gaussian, whose components are statistically independent of the temporal decours of the sources.
  • the source vector s (k) ⁇ is then written as the concatenation of / vectors s ⁇ ⁇ kf. We can therefore write that two components of s (k) are independent if and only if they are associated with two vectors s / k) and sy (k), such that 1 ⁇ j ⁇ f ⁇ J.
  • the mixing matrix A ( ⁇ ) is also written as the concatenation of J matrices A / ⁇ j ), of size (NxPJ), each corresponding to a group of dependent sources.
  • the notations ⁇ and ⁇ j will be omitted.
  • each column vector a ( ⁇ p ) of the matrix A ( ⁇ ) can be written as the product of a gain matrix G ( ⁇ p ), of size (Nx3) , and a vector of quasi-linear parameters (or nuisance) ⁇ p related to the orientation of the tenth source:
  • Vp, 1 ⁇ p ⁇ P, a ( ⁇ p ) G (p p ) ⁇ p (equation 2)
  • ⁇ p the vector of the non-linear and quasi-linear parameters of the /? - th source, where p p is the vector of the nonlinear position parameters and ⁇ ⁇ is the vector of the quasi-linear orientation parameters. It is important to emphasize here that since an application can be based on the same model as that defined by equations (1) and (2), the methods and devices according to the invention are applicable.
  • Each layer 30, 31, 32 represents a different fabric, of conductivity ⁇ j (1 ⁇ j ⁇ 3) assumed to be homogeneous and
  • the p-th dipole is characterized by its position p p , its orientation ⁇ p and its time course s p (k) at time k.
  • F is a so-called switch matrix (according to JC MOSHER, RM LEAHY, and PS 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, May 8-12, 1995, pp. 2943-2946) allowing, in the case where the potentials are recorded by N + 1 sensors, to subtract the value collected by one of the sensors used as reference, to all the values from the N other sensors, and thus to obtain N differences of potentials.
  • the potential created by a dipole in a medium with three homogeneous and isotropic spheres can be approximated by the sum of the potentials created by three dipoles, each dipole being placed in a single homogeneous and isotropic sphere (see Figure 4).
  • the three "unique" spheres are identical to the outer sphere (scalp) of the three-layer model (same radius, same conductivity).
  • the three dipoles are oriented in the same way as the "original" dipole, their positions and their time course are proportional to those of the original dipole (proportionality coefficients ⁇ / for the time decours, ⁇ , - for the positions).
  • This statistical magnitude is calculated from 2q components of the vector x (k), considering q conjugated terms xtfk) and q unconjugated terms.
  • Hermitian C 2 ⁇ x of size (A ⁇ xV), called statistical matrix of order 2q.
  • the matrix C 2 1 ⁇ x can be written, for each 0 ⁇ £ ⁇ q 0 , in the following form: (equation 13) where C, ' (/ iX of size (AZ 1 W), C 2 ⁇ 5 of size (PW) and C 2 ⁇ of size (NxN) are respectively the statistical matrices of order 2q of the vectors x (k), s (k) and v (k).
  • RECTIFIED SHEET (RULE 91) ISA / EP ⁇ (.) is the Kronecker symbol and A ®q is the high Kronecker A matrix where ® represents the Kronecker product operator.
  • the matrix Q x (M) has a structure equivalent to that of the covariance matrix C ⁇ x , which gives UFO-MUSIC authors the possibility of using it directly in MUSIC and RapMUSIC without having to modify these two algorithms. .
  • equation (13) can be rewritten as follows:
  • Jl 1115 is the real matrix of size (R (J, q, €) xR (J, q, €)) of the eigenvalues not
  • E 2q v is the size matrix (iV 9 x (WR (J 5 ⁇ , /)) of the orthonormalized eigenvectors associated with the zero eigenvalues of C 1 1 ⁇ x •
  • C 2q ⁇ x is Hermitian
  • each column vector of E [ qs is orthogonal to each column vector of / J 2 Vv therefore, each column vector of ⁇ (J, q, ⁇ ) is orthogonal to each column vector of E 2 ' qv .
  • the J % criterion can be (equation 17) .
  • the minimization of criterion (17) can be partially achieved by minimizing the following criterion with respect to p: (equation 19) where vpm ⁇ B ⁇ designates the minimal eigenvalue of the matrix B.
  • vpm ⁇ B ⁇ designates the minimal eigenvalue of the matrix B.
  • detfZi ⁇ denotes the determinant of matrix B.
  • RapMUSIC RECTIFIED SHEET
  • RULE 91 RECTIFIED SHEET
  • RULE 91 RECTIFIED SHEET
  • RULE 91 RECTIFIED SHEET
  • RULE 91 RECTIFIED SHEET
  • the RapMUSIC method works with the signal space and therefore applies its orthogonal projector to the signal space. Our process applies our 2nd order projector to the C 2 ⁇ matrix and then calculates the resulting noise space. Our choice to apply the projector before calculating the noise space, aims to further increase the estimation resolution of the parameters of interest.
  • RapMUSIC we propose here a recursive calculation of the projector, which is less expensive in terms of calculation complexity.
  • the latter is obtained by taking the eigenvector associated with the minimum eigenvalue of the matrix
  • (1) can be extracted from ⁇ f (p ⁇ (X) ) •
  • M N q'2 vectors b ⁇ x) ⁇ m) (l ⁇ m ⁇ M) of size (N 2 Xl) (such
  • (1) is estimated by the eigenvector common to all matrices of ⁇
  • ⁇ 1) can be reconstructed only
  • this indeterminacy is not a problem in itself, since, with regard to our current dipole localization problem, the vector ⁇
  • the joint diagonalization can be carried out in our case by means of the method JAD (for "Joint Approximate Diagonalization” in English) proposed by J.-F. CARDOSO and A. SOULOUMIAC in their article entitled “Jacobi angles for simultaneous
  • Another less expensive way, in terms of time and computing resources, is to remove the contribution of the first source estimated not from the observations but rather from the statistical matrix C 2 1 ⁇ x itself.
  • E 2 e qv is increased, which reflects an increase in the size of the noise space, space generated by the column vectors of E l ⁇ qv .
  • the p-th stage of the recurrence of the method, or deflation step, according to the invention consists in minimizing the criterion / 4 of the equation (20) by replacing G q e (p) by ⁇ $ p ) ... ⁇ % 2) ⁇ $ x) G q t (p) and where ⁇ J iK is constructed from the eigenvectors associated with the zero eigenvalues of the matrix ( ⁇
  • ⁇ and ⁇ ⁇ (p) are defined by: and
  • Step 1 Choose the appropriate statistical order 2q according to the number P of sources that one wishes to treat. In practice, q is the minimum value that ensures the processing of all sources potentially present in the multidimensional environment studied.
  • Step 2 Estimate the cumulants C ⁇ ⁇ ] J ' ' ⁇ of order 2q from K samples x (k) and choose the adequate matrix storage for which the estimated statistical matrix of order 2q will be denoted by C ⁇ x (Remember that this statistical matrix is of size (N 9 XN *)).
  • Step 3 (53) Calculate the eigenvalues of the Hermitian matrix C ⁇ x and extract an estimate Ej ⁇ x from the matrix E [ qv .
  • This step may require the estimation of the rank of the matrix C ⁇ x for cases where the number of sources and / or their spatial coherence are / is unknown (s) (hence we will note P the estimate of P).
  • Step 4 (54): Calculate an estimate, ⁇ ⁇ ( ⁇ '' X ) B ⁇ x , of the noise projector ⁇ ⁇ of order 2q.
  • Step 8.2 (58 2 ): Then reassign the variables as follows: - P: ⁇ Pl;
  • Step 9 (59): If P ⁇ N (ie if the source mixture is overdetermined), estimate the P linear parameters s p (k) for each value k by applying to the observation vector x ⁇ k) FAS filter defined by W TcL] -N- ( ⁇ )> OR ⁇ ( ⁇ ) is "a .tau.n at i" i this mixture reconstructed from the estimated ⁇ quasi-linear parameters and non-linear sources using equation (2).

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=l 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

Procédé et dispositif d'identification de paramètres multidimensionnels : application à Ia localisation et Ia 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'électrophysioloeie 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 (électrodes d'enregistrement) disposés de façon standardisée sur le cuir chevelu et
FEUILLE RECTIFIEE (REGLE 91) ISA/EP 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 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 structures cérébrales mises enjeu 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. 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. 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 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 temporelle (<lms), autorisent une étude fine des dynamiques d'activation neuronale. Cependant, l'identification précise des régions cérébrales activées
FEUILLE RECTIFIEE (REGLE 91) ISA/EP 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 et dans l'article de S. MIRON, N. LE BIHAN et J. I. MARS : « Vector-sensor MUSIC for polarized seismic sources localization » paru dans EURASIP
FEUiLLE RECTIFIEE (REGLE 91) ISA/EP Journal on Applied Signal Processing, vol. 10, pp. 74-84, October 2005. L'article précédent présente en détail le contexte spécifique de la sismique d'exploration, et en particulier la formulation du problème direct. 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
FEUILLE RECTIFIEE (REGLE 91) ISA/EP 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 électrode sensitivity - An application of reciprocity », IEEE Transactions On Biomédical 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 Transactions On Biomédical 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 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 Savait, permet de calculer le potentiel électrique et le champ magnétique créés sur les capteurs de surface,
FEUILLE RECTIFIEE (REGLE 91) ISA/EP 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 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 Elément
Method, en anglais), bien décrite dans l'article de A. S. FERGUSON, X. ZHANG, et G. STROINK, "A complète linear discretization for calculating the magnetic field using the boundary élément method," IEEE Transactions On Biomédical Engineering, vol. 41, pp. 455-459, 1994 ;
- la méthode des éléments finis (FEM, pour Finite Elément 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 Différence 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 différences » Phys. Med. Biol., vol. 41, no. 7, pp 1079-1091, 1996.
FEUILLE RECTIFIEE (REGLE 91) ISA/EP 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 Conférence on Acoustics Speech and Signal Processing, vol. 5, Détroit, Michigan, mai 1995, pp. 2943-2946 et : « EEG and MEG: Forward solutions for inverse methods », IEEE Transactions
On Biomédical Engineering, vol. 46, no. 3, pp. 245-259, Mars 1999, ont proposé un mode de calcul commun basé sur i) une formulation matricielle spatiotemporelle 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, et plus généralement les paramètres non linéaires et quasi-linéaires des sources d'intérêt.
L'une des plus célèbres de ces méthodes est la méthode MUSIC (MUltiple
Signal Classification en anglais, ou classification de signaux multiples, en français) exploitant les statistiques d'ordre 2 et la matrice de covariance des observations, décrite par R. O. SCHMIDT dans ses publications : « A signal subspace approach to multiple emitter location and spectral estimation », thèse de
FEUILLE RECTIFIEE (REGLE 91) ISA/EP doctorat de l'UIniversité 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 Porthogonalité 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éeselle ne tire pas avantage de la possible factorisation du problème direct. En outre, bien qu'en théorie la méthode MUSIC permet d'estimer à la fois les paramètres non linéaires et les paramètres quasi-linéaires, ceci est difficilement réalisable en contexte opérationnel de par la complexité de calcul que cela entraine
II faut 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 réaliser en contexte pratique. En effet, E. FERRARA propose une nouvelle version de MUSIC qui exploite la factorisation de la formulation matricielle du problème direct, et dissocie ainsi l'estimation des paramètres non-linéaires (c'est- à-dire les positions des dipôles de courant en électrophysiologie humaine) de celle des paramètres quasi-linéaires (c'est-à-dire les orientations des dits dipôles) des sources, réduisant de ce fait considérablement le coût de calcul de l'algorithme. L'algorithme de E. FERRARA rend alors réalisable en contexte pratique l'estimation simultanée des paramètres non linéaires et quasi-linéaires.
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 »,
/EEE Transactions On Signal Processing, vol. 47, no. 2, pp. 332-340, Février
FEUILLE RECTIFIEE (REGLE 91) ISA/EP 1999. Cette dernière permet notamment de faciliter la recherche de plusieurs optima locaux, ou en d'autres termes d'améliorer l'estimation des paramètres d'intérêt dans la métrique MUSIC, en particulier lorsque la dimension de l'espace source augmente ou bien lorsque les paramètres non linéaires sont très voisins d'une source à l' autre.
Contrairement aux autres approches séquentielles, la méthode RapMUSIC s'appuie sur les travaux de E. FERRARA, offrant ainsi la possibilité d'exploiter une éventuelle factorisation de la formulation matricielle du problème direct en présence d'une antenne à diversité de polarisation, par exemple, et plus généralement lorsque les paramètres quasi-linéaires sont inconnus et doivent êtres estimés.
Alors que les approches précédentes n'exploitent que les statistiques d'ordre 2, et donc la matrice de covariance des observations, B. PORAT et B. FRIEDLANDER décident de proposer une version de MUSIC exploitant cette fois les statistiques d'ordre 4, et donc la matrice de quadricovariance des observations, tel que décrit dans leur article intitulé « Direction finding algorithme 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 et est différente algorithmiquement de la version originale de MUSIC non seulement de par l'utilisation des statistiques d'ordre supérieur, mais aussi de par le fait que la structure algébrique de la matrice de quadricovariance est différente de celle de la matrice de covariance, et demande alors aux auteurs de modifier la version originale de MUSIC.
Plus récemment, SATOSHI NIIJIMA et SHOOGO UENO ont proposé une méthode, baptisée UFO-MUSIC, exploitant à la fois les statistiques d'ordre 4 et une éventuelle factorisation de la formulation matricielle du problème direct, tel
FEUILLE RECTIFIEE (REGLE 91) ISA/EP que décrit dans l'article « Universal Fourth Order Music Method: Incorporation of ICA into MEG Inverse Solution », Third international conférence on Independent Component Analysis and Blind Signal Séparation, 9-12 décembre 2001, San Diego, California, USA. Cependant, si les auteurs de cette méthode disent exploiter les statistiques d'ordre 4, ils les exploitent sous une fome différente de celle exploitée par B. PORAT et B. FRJEDLANDER. En effet, ils exploitent une ou plusieurs matrices dites de quadricovariance contractée dont la structure algébrique est différente de celle de la matrice de quadricovariance de B. PORAT et B. FRIEDLANDER. Mais est par contre équivalente à celle d'une matrice de covariance, ce qui leur permet d'utiliser l'algorithme MUSIC de E. FERRARA ou même la méthode RapMUSIC sans modification de ces derniers lorsqu'une seule matrice de quadricovariance contractée est utilisée. Lorsque plusieurs matrices sont utilisées, les auteurs proposent d'estimer l'espace signal et l'espace bruit par diagonalisation conjointe des différentes matrices de quadricovariance contractée employées. Les algorithmes MUSIC et RapMUSIC peuvent alors à nouveau être utilisés sans qu'il y ait besoin de les modifier.
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 », /EEE 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, quasi 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 séparation of independent narrow-band sources: Concept and
FEUILLE RECTIFIEE (REGLE 91) ISA/EP 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 exclusivement la localisation et la reconstruction d'activités électriques intracérébrales, d'autres algorithmes ont vu le jour. Le lecteur pourra se référer à l'article de C. M. MICHEL, M. M. MURRAY, G. LAΝTZ, 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 ou lorsqu'un certain niveau de résolution est requis.
FEUILLE RECTIFIEE (REGLE 91) ISA/EP 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 - I 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. Par conséquent, une implémentation de cet algorithme n'est pas réalisable en contexte opérationnel lorsque les paramètres quasi-linéaires sont inconnus et doivent êtres estimés. De plus, cette méthode souffre des problèmes justifiant l'emploi des approches séquentielles. Remarquons qu'une évolution de cette méthode vers un algorithme séquentiel, qui plus est, exploitant une éventuelle factorisation du problème direct n'est en aucun cas trivial car comme nous l'avons écrit la matrice de quadricovariance de B. PORAT a une structure algébrique différente de celle de la matrice de covarariance. En outre cette évolution n'a jamais été proposée ni même suggérée dans l'art antérieur. Observant à présent la méthode UFO-MUSIC proposée par SATOSHI NIIJIMA et SHOOGO UENO, si l'idée d'exploiter conjointement l'information contenue dans les différentes matrices de quadricovariance contractée est intéressante en soi, la mise en œuvre proposée par les auteurs n'est pas justifiée. En effet, l'algorithme de diagonalisation conjointe utilisé cherche une matrice de passage orthonormée commune aux différentes matrices de quadricovariance contractée utilisées, mais l'existence d'une telle solution n'est pas démontrée par les auteurs. Pour cause, sauf cas particuliers non exploitables en contexte opérationnel, rien ne garantit l'existence d'une telle solution.
En outre, l'approche proposée par SATOSHI NIIJIMA et SHOOGO UENO ne permet pas de résoudre de problème inverse mal posé tel que rencontré lorsque le nombre de sources, P, est supérieur au nombre d'observations, N. Ceci est principalement dû au fait que les auteurs préfèrent exploiter conjointement une
FEUILLE RECTIFIEE (REGLE 91) ISA/EP ou plusieurs combinaisons linéaires de tranches matricielles du tenseur cumulant d'ordre 4 plutôt que le tenseur lui-même.
De plus, SATOSHI NIUIMA et SHOOGO UENO supposent dans la mise en œuvre de leur méthode que les sources sont spatialement statistiquement indépendantes à l'ordre 4, ce qui en contexte opérationnel est loin d'être toujours vérifié. Par exemple, en épilepsie, deux zones épileptiques peuvent avoir des activités électriques fortement corrélées et pourtant se trouver distantes de plusieurs centimètres.
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 « 2g-RapMUSIC »
(q≥2), qui est la contraction de « méthodes 2d-D-MUSIC », 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
FEUILLE RECTIFIEE (REGLE 91) ISA/EP 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 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
FEUILLE RECTIFIEE (REGLE 91) ISA/EP 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 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), et/ou, contrairement à la méthode UFO-MUSIC, la structure algébrique du tenseur cumulant d'ordre 2q (q ≥ 2) dans sa globalité et non celle de tranches matricielles du dit tenseur, afin de pouvoir simultanément i) résoudre un problème mal posé en ce 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 de l'invention vise à offrir une telle technique capable d'exploiter une éventuelle factorisation de la formulation matricielle du problème direct afin de dissocier l'estimation des paramètres non linéaires des sources de l'estimation des paramètres quasi-linéaires, et de réduire ainsi sensiblement les coûts de calcul lorsque les paramètres quasi-linéaires sont inconnus. Cet objectif vise donc à garantir la faisabilité du procédé en contexte opératiσnel lorsque les paramètres quasi-linéaires des sources sont inconnus. Un objectif supplémentaire est de fournir une telle technique qui permette d'exploiter le concept de déflation étendu et formalisé à 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, notamment dans l'espace.
FEUILLE RECTIFIEE (REGLE 91) ISA/EP Encore un objectif de l'invention vise à fournir une telle technique n'exploitant aucune autre information a priori sur les paramètres linéaires des sources.
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 œuvre, y compris dans des dispositifs, par exemple du type dispositif d'électroencéphalographie, de magnétoencéphalographie, d'étude ou d'exploration 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 iV≥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.
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≥Z) des observations, cette matrice étant de taille (N^dV) et, pour q≥2, de structure algébrique différente de celle de la matrice de covariance ;
FEUJLLE RECTIFIEE (REGLE 91) ISA/EP - 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 œuvre une étape de localisation et de reconstruction de l'activité électrique générée par la pluralité de P≥i 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(Θ)s(k) + v(k) où : - s(!c) est un vecteur, de taille (Px 1), représentatif des paramètres linéaires correspondant aux décours temporels desdites P sources d'intérêt, non
FEUILLE RECTIFIEE (REGLE 91) ISA/EP gaussiennes et potentiellement (mais pas totalement) corrélées suivant ledit au moins un premier paramètre multidimensionnel ;
- A(Θ) est une matrice de mélange instantané, de taille (NxP), où <9 = {0, θp} 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(Θ) se décompose sous la forme : a(θp) ≈G(pp)Φp avec pp et φp représentant respectivement les paramètres non linéaires d'une part et quasi-linéaires d'autre part 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 (NxI), du bruit additif, indépendant des sources d'intérêt.
Préférentiellement, le procédé selon l'invention comporte en outre au moins une étape i) d'estimation des cumulants Cî"*1"; ";» 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 (WxN*1), sera notée C^x .
De façon avantageuse, le procédé selon l'invention comporte également au moins une étape d'estimation du rang de la matrice C^x , 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 C^x et une étape de construction d'une fonction de coût, dite pseudo-spectre d'ordre 2q ou 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>)H π;,,G,>)} % forme JAp) = — —, z ^ , ou άet{G<(p)HG;(p)}
FEUILLE RECTIFIEE (REGLE 91) ISA/EP - ήj* = (JS2^x ) Ê[ζx est un opérateur matriciel, dit projecteur bruit d'ordre 2q, avec Ê2 e*v la matrice des vecteurs propres orthonormalisés associés aux valeurs propres nulles de ladite matrice C^x ;
- G^ (/7) =G!(p)®î"{ ®G(p)*®f 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 φ. De façon avantageuse, le procédé selon l'invention met aussi en œuvre 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 des P sources d'intérêt.
Préférentiellement, la /?-ième étape (l≤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 p^p) ;
- calcul d'un vecteur Φq °μ (Pξ(P)) en prenant le vecteur propre correspondant à la valeur propre minimale de la matrice
[G:°- (AW)V UW)]V1 (%))H KWÏ (AM) •»
- extraction d'un vecteur ^ représentatif d'une estimée du vecteur de paramètres quasi-linéaires ^(p) , à partir du vecteur φ*°>° (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) construction d'un vecteur α (*%;>) ) = G(%>))4(/>) ' Puis ") calcul d'une matrice ij^ de taille (Λ^x/V9) tenant compte d'un remplacement dudit vecteur a[θξipλ par ledit vecteur o(θξ, λ et iii) réaffectation des variables suivant les fonctions suivantes :
• P := P~ \ ;
FEUILLE RECTIFIEE (REGLE 91) ISA/EP De façon préférentielle, l'étape d'extraction d'un vecteur 4^ comprend les sous-étapes suivantes :
- extraction de M = N*'2 vecteurs bξ qfy (m) de taille (N2X 1) ;
- transformation des vecteurs en M matrices Bξ qfτ (m) de taille (NxN) ; - calcul d'un vecteur propre commun aux M matrices de l'ensemble
Δ*^' 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 œuvre une étape de construction d'un filtre FAS (pour « filtrage Alterné Séquentiel ») défini par la formule par W = [(?"__ J 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 processeur adapté à mettre en œuvre 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≥\ 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é. 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 magnétoencéphalographe dans le domaine neurologique.
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 œuvre 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
FEUILLE RECTIFIEE (REGLE 91) ISA/EP 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, 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≥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≥l 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 ;
FEUILLE RECTIFIEE (REGLE 91) ISA/EP - 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
FEUILLE RECTIFIEE (REGLE 91) ISA/EP 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. Pour pouvoir appliquer ce mode de réalisation à la MEG ou bien à l'exploration sismique, il est nécessaire de prendre connaissance de la formulation du problème direct propre à ces deux domaines d'apllications. De tels renseignements sont accessibles dans l'ouvrage de BIN HE intitulé « Modeling and Imaging of Bioelectrical Activity. Principles and Applications » paru en 2004 aux éditions KLUMER ACADEMIC/PLENUM PUBLISHERS, NFEW YORK en ce qui concerne la MEG, et dans l'article de S. MIRON, N. LE BIHAN et J. I. MARS : « Vector-sensor MUSIC for polarized seismic sources localization » paru dans EURASIP Journal on Applied Signal Processing, vol. 10, pp. 74-84,
October 2005 pour ce qui est de l'exploration sismique. 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 iV+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(θ)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(Θ) est la matrice de mélange instantané de taille (NxP), où O=[U1,...,θP) désigne l'ensemble des paramètres quasi-linéaires et non linéaires associés aux P sources, Le., 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.
FEUILLE RECTIFIEE (REGLE 91) ISA/EP 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. On notera Pj le nombre de sources du y-ème groupe, et P = £"!_, Pj . Notons que le cas où J=P correspond au cas où les P sources sont statistiquement indépendantes alors que le cas où J=J correspond à celui où elles sont toutes corrélées.
Le vecteur des sources s(k)Υ s'écrit alors comme la concaténation de / vecteurs s}{kf. On peut donc écrire que deux composantes de s(k) sont indépendantes si et seulement si elles sont associées à deux vecteurs s/k) et sy(k) , tels que 1 ≤j ≠f ≤ J.
La matrice de mélange A(Θ) s'écrit également comme la concaténation de J matrices A/Θj), de taille (NxPJ), correspondant chacune à un groupe de sources dépendantes. Dans la suite, et par souci de simplicité, les notations Θ et Θj seront omises.
Dans le cadre de la localisation de sources en EEG, chaque vecteur colonne a(θp) de la matrice A(Θ) peut s'écrire comme le produit d'une matrice de gain G(θp), de taille (Nx3), et d'un vecteur de paramètres quasi-linéaires (ou de nuisance) Φp lié à l'orientation de lap-ième source :
Vp, 1 < p ≤ P, a(θp) =G(pp)Φp (équation 2)
On notera θp
Figure imgf000026_0001
le vecteur des paramètres non linéaires et quasi- linéaires de la /?-ième source, où pp est le vecteur des paramètres non-linéaires de position et Φμ le vecteur des paramètres quasi-linéaires d'orientation. II 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é σj ( 1 < j ≤ 3 ) supposée homogène et
FEUILLE RECTIFIEE (REGLE 91) ISA/EP 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 Φp et son décours temporel sp(k) a 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 multiple-shell spherical head models »
Electroencephαlogrαphy αnd Clinicαl Neurophysiology, vol. 90, no. 1, pp. 58-64, January 1994, et illustrée dans la figure 4 :
G(Pp) (équation 3)
Figure imgf000027_0001
où F 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 Conférence on Acoustics Speech and Signal Processing, vol. 5, Détroit, 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é λ/ pour les décours temporels, μ,- pour les positions).
FEUILLE RECTIFIEE (REGLE 91) ISA/EP Le vecteur h(r,p), de taille (Nx3) est donné quant-à-lui par l'expression suivante : h(r,p) = [(c, -c2
Figure imgf000028_0001
(équation 4) où les paramètres c\ et et sont définis par :
5.)
Figure imgf000028_0002
avec σ la conductivité de la sphère extérieure (scalp) et :
-y-Aiw-A+yt-f') (équation 6)
Notons que {λu A2, A3, μ\, μ2, //3} 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 sphères », Physics in Medicine and Biology, vol. 40, no. 3, pp. 335-349, Mars 1995, en minimisant la quantité suivante :
Δ (équation 7)
Figure imgf000028_0003
où le nombre de termes iVmax doit être suffisamment grand pour garantir une bonne précision de calcul (par exemple Nmax = 200), q et Rj ( 1 < y < 3 ) sont,
respectivement, les conductivités et les rayons des trois sphères, et :
Vu, /_ (équation 8)
Figure imgf000028_0004
FEUILLE RECTIFIEE (REGLE 91) ISA/EP Les coefficients ni2\ et W22 sont alors obtenus comme suit :
)
Figure imgf000029_0001
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 Biomédical
Engineering, vol. 46, ηo. 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 Elément Method), telles que représentées en figure
3.b.
6.2 Statistiques d'ordre 2a
On considère les statistiques d'ordre Iq des observations. On définit, pour tout k, le cumulant d'ordre 2q d'un vecteur x(k) de la façon suivante :
Figure imgf000029_0002
(équation 10)
Cette grandeur statistique est calculée à partir de 2q composantes du vecteur x(k), en considérant q termes xtfk) 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 Cuy*' .
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, Spécial Issue on Linear Algebra in Signal and Image Processing, vol. 39 IC, pp.
3-30, November 2004.
FEUILLE RECTIFIEE (REGLE 91) ISA/EP 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 C2^x , de taille (A^x-V), appelée matrice statistique d'ordre 2q.
Plusieurs façons de ranger ces statistiques d'ordre 2q dans la matrice C2^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
#o=l.
Ces qo arrangements sont indexés suivant l'entier € (O≤€≤qo-1) et correspondent respectivement à qo matrices statistiques d'ordre 2q, notées Cf,%x . Le (//,/,' )-ième élément (l ≤ I',I2 ≤ Nq ) de la matrice C2 1^x est alors donné par l'expression suivante :
<., (>/>'2) = Çb;J« (équation 11) où pour chaque nombre € tel que O≤€≤qo-l et pour chaque indice (,- tel que l ≤
(équation 12)
Figure imgf000030_0001
La propriété de multilinéarité des cumulants et des moments permet d'écrire chaque matrice C2 1^x {q ≥ 1 ) sous une forme spéciale.
En effet, sous l'hypothèse d'indépendance des sources et du bruit, la matrice C2 1^x peut s'écrire, pour chaque 0< £ ≤ q0 , sous la forme suivante : (équation 13)
Figure imgf000030_0002
où C,'(/iX de taille (AZ1W), C2^5 de taille (PW) et C2^ de taille (NxN) sont respectivement les matrices statistiques d'ordre 2q des vecteurs x(k), s(k) et v(k).
FEUILLE RECTIFIEE (REGLE 91) ISA/EP δ(.) est le symbole de Kronecker et A®q désigne la matrice A élevée à la puissance de Kronecker où ® représente l'opérateur du produit de Kronecker. Le nombre I est le même que celui défini dans les équations (11) et (12). Il est important de remarquer que pour q≥2 la structure algébrique de la matrice statistique
Figure imgf000031_0001
exploitée dans les algorithmes MUSIC et RapMUSIC. En outre, la C^q x pour q=2 est différente de la matrice de quadricovariance contractée, Qx(M), exploitée dans la méthode UFO-MUSIC. En effet si, pour q=2 et en présence de données à valeurs réelles, la matrice C^q x de taille (N2xN2) admet d'après l'équation (13) la structure suivante :
Figure imgf000031_0002
la matrice de quadricovariance contractée, Qx(M), de taille (NxN), admet quant à elle la structure suivante :
Qx(M) =Me) Qs(B)A(Oj" avec :
IQs(B)IiJ = ∑Λ*=-..P CUm(S1-, Sy, Sn Ss,) Brs et :
Figure imgf000031_0003
où α(θr) désigne Ie r-ième vecteur colonne de la matrice de transfert A(θ) et où la matrice Q8(B) est de taille (PxP). Ainsi la matrice Qx(M) a une structure équivalente à celle de la matrice de covariance C^x ce qui donne aux auteurs de UFO-MUSIC la possibilité de l'utiliser directement dans MUSIC et RapMUSIC sans avoir besoin de modifier ces deux algorithmes.
D'autre part, on note que l'équation (13) peut être réécrite de la façon suivante :
<.* = ∑[Af~e ® A™]tfφSj[Af>-< ®Aγ]H +δ{q-
Figure imgf000031_0004
(équation 14), où la matrice Cf,/ S , de taille ( PJ x PJ), est la matrice statistique d'ordre 2q du vecteur sβ).
U 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)
FEUILLE RECTIFIEE (REGLE 91) ISA/EP à 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, Spécial 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 « 2^-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) V l≤ y ≤ Λ Pj < N ;
A2) La matrice Af"'' ® Af est de rang plein Pf ;
A3) R(J,q,i) = ∑'HP? < N-> ;
A4) La matrice A[J^J) = [Af"'' ®Af,...,Af-e 0Af] est de rang plein
R{J,q,Q. En particulier, pour P sources statistiquement dépendantes, c'est-à-dire pour lesquelles J =1, les hypothèses (A1)-(A2) se réduisent donc à : A'1) P < N ;
A'2) la matrice A(IqJ) ≈ A®"'1 ® Aw est de rang plein P". Alors que pour P sources statistiquement indépendantes, c'est-à-dire pour lesquelles J=P, les hypothèses (A1)-(A2) se réduisent alors à:
FEUILLE RECTIFIEE (REGLE 91) ISA/EP A"l) P < N" ;
A"2) La matrice Λ(P,q,t) ≈ A0<l'ι0At0t est de rang plein égal à P ; où 0 et A09 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 (q ≥ 2 ), ainsi que l'approche dénommée ici par les inventeurs « 2^-RapMUSIC ». 6.3.1 - Pseudo-spectre d'ordre 2a (a≥2 )
A partir de l'équation (13), il est désormais possible de calculer la décomposition en valeurs propres de la matrice hermitienne C2q x :
C=O^ K,] ^S I [K, *.'„]" (équation 15)
où Jl1115 est la matrice réelle de taille (R(J,q,€)xR(J,q,€)) des valeurs propres non
nulles de
Figure imgf000033_0001
des vecteurs propres
orthonormalisés associés.
D'autre part, E2q v est la matrice de taille (iV9 x( W-R(J5^,/)) des vecteurs propres orthonormalisés associés aux valeurs propres nulles de C1 1^x • Comme la matrice C2q<x est hermitienne, chaque vecteur colonne de E[q s est orthogonal à chaque vecteur colonne de /J2Vv
Figure imgf000033_0002
, par conséquent, chaque vecteur colonne de Λ(J,q,έ) est orthogonal à chaque vecteur colonne de E2'q v .
Ainsi, en notant θp s le vecteur des paramètres de localisation (position et orientation) de la /?-ième (l≤p≤P) source, appartenant qui plus est au jf-ième groupe, et a (θp J f~' ® P. ) "-ième colonne de la matrice
Figure imgf000033_0003
(l≤p≤P et l ≤ j ≤ J ) sont tous orthogonaux aux vecteurs colonne de E2qΛ, . Ce résultat nous donne la possibilité de construire une fonction de coût dont la minimisation
FEUILLE RECTIFIEE (REGLE 91) ISA/EP permet d'estimer les P vecteurs θj, . Cette fonction de coût, baptisée pseudospectre d'ordre 2q (ou pseudo-poly spectre), est définie par : )•*"' ®α(0)*β']H πlv[a(θfq-1 ®a(θrι] (équation 16) où Uq ι v
Figure imgf000034_0001
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 {θ) = J1 (θ)/L(θf~c ®a(θ)^ au lieu de celui défini par l'équation 16, afin de rendre le pseudo-poly spectre constant par rapport à θ en l'absence de sources.
Puis, d'après l'équation (2) et les propriétés du produit .de Kronecker, le critère J% peut être (équation 17)
Figure imgf000034_0002
Figure imgf000034_0003
.
Minimiser le critère J2 de l'équation (17) par rapport à θ revient à i) minimiser la valeur propre minimale, notée ^1 {p) , de la matrice Gq' (pf uq vGq ( (p) dans la métrique
Figure imgf000034_0004
par rapport à p tel que papt = arg min {/Ij (/>)} , et ii) déduire de popt le vecteur Φq* solution de la p minimisation du critère T2 par rapport à θ et notéΦq ι [pBpl) , en prenant le vecteur propre associé à la valeur propre
Figure imgf000034_0005
.
En effet, P&,{p<>p<) et φé{p»pι) vérifient :
Figure imgf000034_0006
(é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 : (équation 19)
Figure imgf000034_0007
où vpm{B} désigne la valeur propre minimale de la matrice B. Nous passons ainsi d'une optimisation dans un sous-espace vectoriel de R6 , à une optimisation dans un sous-espace vectoriel de E3 , 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 :
FEUILLE RECTIFIEE (REGLE 91) ISA/EP (équation 20)
Figure imgf000035_0001
où detfZï} 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 /3 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 Φ* 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 J$ 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 Φ? pour toutes valeurs de q et de €.
6.3.2 - La méthode 2q-RapMUSlC ( a > 2 )
On décrit à présent la méthode dénommée « 2#-RapMUSïC » ( 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 et formalisé à l'ordre 2q et prenant en compte la présence de sources potentiellement corrélées.
Il est important de remarquer que le principe de déflation présenté ici ne peut être interprété comme celui de la méthode RapMUSIC pour laquelle la matrice de covariance aurait été remplacée par notre matrice statistique C^ x d'ordre 2q, notamment parce que la différence de structure algébrique entre les deux matrices ne permet pas d'utiliser le projecteur orthogonal défini dans RapMUSIC. De plus, comme nous le verrons dans cette section, la définition du
FEUILLE RECTIFIEE (REGLE 91) ISA/EP projecteur orthogonal d'ordre 2q n'a rien de trivial, notamment lorsque les sources d'intérêt sont spatialement corrélées. En outre, la méthode RapMUSIC travaille avec l'espace signal et applique donc son projecteur orthogonal à l'espace signal. Notre procédé applique pour sa part notre projecteur d'ordre 2q à la matrice C2^ pour ensuite calculer l'espace bruit résultant. Notre choix d'appliquer le projecteur avant calcul de l'espace bruit, a pour but d'accroître davantage la résolution d'estimation des paramètres d'intérêt. Enfin, contrairement à RapMUSIC, nous proposons ici un calcul récursif du projecteur, ce qui est moins coûteux en termes de complexité de calcul. 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. 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 J1-J4 non pas P minima globaux correspondant respectivement aux P sources, mais plutôt P-I 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 ξ{.) 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'œil 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(Θ) correspondantes sont rangées ne change pas l'expression de x(k).
FEUILLE RECTIFIEE (REGLE 91) ISA/EP Le vecteur φ*{pξ{^ , est quant-à-lui obtenu comme le vecteur normalisé qui, multiplié à gauche par la matrice G q {Pφ)) donne le vecteur a(θm) ®a(θφ)) .
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
[G:(^))HG:(%)J1G:(%))HΠ:.G:(%))-
Puis, le vecteur <|(1) peut être extrait de φf (pξ(X)) • II faut pour cela extraire en premier lieu les M =Nq'2 vecteurs b^x) {m) (l ≤ m ≤ M ) de taille (N2Xl) (tels
que 4e
Figure imgf000037_0001
. - .6|{'f) (M )τ] ), puis les convertir en M matrices
Figure imgf000037_0002
(Ht) est faite des N éléments consécutifs de ^ (m) à partir du [N(«-l)+l]-ième), et enfin diagonaliser conjointement l'ensemble Δ1^,) de matrices définies par {*$) (m) Bfa (m)H , fi*J (m)τ Bfa (m)* , tel que l≤m≤M, si €≈O, par et sinon par
Figure imgf000037_0003
Plus précisément, le vecteur ^|(1) est estimé par le vecteur propre commun à toutes les matrices de Δ|'( c,), et associé à la plus forte valeur propre, est, à un facteur d'échelle près, égal à (|(l) . Nous parlons d'estimation car le vecteur (|{1) 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 <|(1) est un vecteur d'orientation. De ce fait, la connaissance du vecteur directeur, de norme unité, de (|(1) 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
FEUILLE RECTIFIEE (REGLE 91) ISA/EP 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 sensé with application in blind source séparation », IEEE Transactions On Signal Processing, vol. 50, no. 7, pp. 1545-1553, JuIy 2002.
Il est donc facile d'obtenir le vecteur de positions pξ{x), de même que le vecteur d'orientation correspondant. Il est cependant plus difficile de trouver avec précision les P-I minima locaux de J 4 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. Notons que par « mise en œuvre », on entend notamment une formalisation à l'ordre 2q du concept de déflation.
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, simple et triviale, n'en est pas moins très coûteuse car elle nécessite d'estimer P fois la matrice statistique des observations C2q_x , en contexte opérationnel.
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 l'espace signal de la matrice statistique C2q_x • Cette solution découle naturellement de l'approche employée dans RapMUSIC et impose de travailler avec l'espace signal plutôt qu'avec l'espace bruit. Toutefois, cette solution ne laisse pas la possibilité de travailler avec l'espace bruit.
Une autre 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 C2 1^x elle-même.
FEUILLE RECTIFIEE (REGLE 91) ISA/EP Cette approche nécessite d'étudier précisément la structure algébrique de la matrice C2'9 X , de façon à annuler les colonnes de la matrice A®*"' ® A*®' impliquant le vecteur a[0ξ{ι)) = G( pξ{x) J (|(1) à 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 (N&xl), (Ncxl), (Λ/rfXl), le produit de Kronecker vérifie alors b ®c ®d =(lNb ®c ®lNj )(b ®d) .
En effet, cette propriété permet de montrer la proposition suivante pour deux matrices Ψ^,) et Ξ^(l) , de taille (Λ^xiV^) définies par, pour tout entier me {θ,l,...,q-\] :
W«.-l _ T f(l) = IΛT
Figure imgf000039_0001
-jt.m T m t. m (mt.m \ mt.m /m'.'" \
Encore une autre manière de procéder, également moins coûteuse, différente de l'approche utilisée dans RapMUSIC à l'ordre 2, et préférée des inventeurs de la présente demande, consiste tout d'abord à retirer de C1 1^x Ia contribution de la première source estimée, puis à calculer l'espace bruit correspondant à la nouvelle matrice statistique résultante. Ce procédé permet à chaque déflation d'accroître la dimension de l'espace bruit et donc d'améliorer l'estimation des paramètres de la source recherchée.
Cette approche nécessite d'étudier précisément la structure algébrique de la matrice C^x , notamment au travers de l'équation (13). En effet, l'équation (13) nous apprend que nous serons capables de retirer retirer de C^q x la contribution de la première source estimée à condition de pouvoir annuler les vecteurs colonnes de la matrice A®"~f ® A'Θt impliquant le vecteur Λ (^(1) ) = G(/^(1) W{1) , autrement dit les vecteurs colonnes de A®"~e ®A*®e de la forme
Figure imgf000039_0002
ou encore 6®α(^(i>) ®c Ceci n'est en rien trivial et ne peut être déduit de la procédure de déflation de RapMUSIC.
FEUILLE RECTIFIEE (REGLE 91) ISA/EP Ainsi il nous faut construire la matrice 3^(1) de taille (NxN) définie à
partir du vecteur colonne α(6|(i)) de la manière suivante :
Figure imgf000040_0001
Alors, multiplier la matrice Λ®?~' ®A'®' de taille (N9X/3*) à gauche par Σ|J,J
Figure imgf000040_0002
permet d'annuler tous les vecteurs colonnes de
A®q-t Θ A*®ι jmpiiquant a[θm) .
Le vecteur de position, pξ(i) , de la df(2)-ième source est alors obtenu en minimisant le critère J4 de l'équation (20), dans laquelle Gq (p) est désormais remplacé par Σ^q Gq e (p) et où π^ ne sera plus obtenu à partir de la diagonalisation de C^x , mais plutôt à partir de la diagonalisation de
On remarque alors que le rang de la matrice Σ.^C^ (∑|^ ) est à présent strictement plus petit que R(J,q,€). En effet, nous avons diminué le rang de C2 e 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
E2 e q v est augmenté, ce qui traduit une augmentation de la dimension de l'espace bruit, espace engendré par les vecteurs colonnes de El ι q v .
En rappelant que la matrice E{q V intervient directement dans le calcul du projecteur bruit π£ „ = {E[q v ) (E^) , et donc dans le calcul du critère J 4 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 vecteur d'orientation de la £(2)-ième source, ^(2) , est extrait en suivant la procédure détaillée précédemment pour ^(1) et le vecteur directeur de la £(2)-ième source, a f θ^ 1 = G ( /^2) ) 4(2) Peut a*ors facilement être reconstruit.
FEUILLE RECTIFIEE (REGLE 91) ISA/EP II est ensuite procédé par récurrence jusqu'à l'estimation des P vecteurs de paramètres θp
Figure imgf000041_0001
Φp τ~] correspondant aux P sources d'intérêt.
On précise également que la p-ième étape de la récurrence du procédé, ou étape de déflation, selon l'invention consiste à minimiser le critère /4 de l'équation (20) en remplaçant Gq e (p) par ∑$p) ...Σ%2)Σ$x)Gq t(p) et où πJiK est construit à partir des vecteurs propres associés aux valeurs propres nulles de la matrice
Figure imgf000041_0002
(∑|('p)...Σ»j'2)Σ|(())H . Selon l'encore autre manière de procéder les matrices ∑|^ et Ξ^(p) sont définies par :
Figure imgf000041_0003
et :
Contrairement à l'opérateur de déflation proposé à l'équation (13) de l'article décrivant RapMUSIC, la construction de notre opérateur est récursive et ne nécessite ainsi pas d'inversion de matrice : le coût de calcul s'en trouve alors fortement réduit.
6.3.3 - IdenMabïlité
On peut montrer que le problème du traitement à l'ordre Iq 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 €, C^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 A/lq sont différents. Il est important de noter que le nombre Λ/£ dépend implicitement de la structure algébrique du vecteur a(θ) = G(p)Φ , 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 t fixées, le nombre Λ/f 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 «
Figure imgf000041_0004
où d est la dimension du vecteur Φ, autrement dit le nombre de paramètres quasi- linéaires propres à une source.
FEUILLE RECTIFIEE (REGLE 91) ISA/EP 6.4 Résumé des grandes étapes successives du procédé selon l'invention
Relativement à la figure 5, on résume ci-dessous les différentes étapes du procédé selon l'invention en supposant que K échantillons temporels d'observations de surface, x(k) (l ≤ 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ζι ]J''^ 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 C^x (on rappelle que cette matrice statistique est de taille (N9XN*)).
Etape 3 (53) : Calculer les valeurs propres de la matrice hermitienne C^x et extraire une estimée Ej^x , de la matrice E[q v . Cette étape peut nécessiter l'estimation du rang de la matrice C^x 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, π^ = (Ê^''X) B^x , du projecteur bruit π^ d'ordre 2q.
Etape 5 (55) : Calculer une estimée, J4 , du critère J4 de l'équation 20 en utilisant la matrice π^ , pour une grille de positions appropriée selon les termes précités, et chercher son minimum global, noté /L1, .
Etape 6 (56) : Calculer le vecteur φ^m [pξ<Λ en prenant le vecteur propre correspondant à la valeur propre minimale de la matrice
[Gr' M" Qrt (pm )] V' (p mf iW" (Pm) •
Etape 7 (57) : Extraire le vecteur <|{Ij , estimée du vecteur quasi-linéaire
<|(1) , à partir du vecteur $"* (/â(1)) • Pour cela, extraire tout d'abord les M = ISf'2 vecteurs
Figure imgf000042_0001
de taille (N2X 1), puis les transformer alors en M matrices B^'" (m) de taille (NxN), et enfin calculer le vecteur propre commun aux M
FEUILLE RECTIFIEE (REGLE 91) ISA/EP matrices de l'ensemble λt °,p' 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 (58i) : 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
Figure imgf000043_0001
de la manière décrite dans la section précédente en remplaçant le vecteur a{θξ,Λ par
Etape 8.2 (582) : Puis, réaffecter les variables de la manière suivante : - P :≈ P-l ;
Figure imgf000043_0002
et retourner à l'étape 3 (53), Sinon aller à l'étape 9 (59). 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 W =TcL] -A(^) > OU ^ (^) est "a τnati"ice de mélange reconstruite à partir de l'estimation Θ des paramètres quasi-linéaires et non linéaires des sources en utilisant l'équation (2).
FEUILLE RECTIFIEE (REGLE 91) ISA/EP

Claims

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≥\ 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≥\ 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) de taille (Λ^x/V7) ;
- 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 œuvre une étape de localisation et de reconstruction de l'activité électrique générée par ladite pluralité de P≥\ 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ètres
FEUILLE RECTIFIEE (REGLE 91) ISA/EP 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 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(Θ)s(k) + v{k) où :
- s(k) est un vecteur, de taille (PxI), 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(Θ) est une matrice de mélange instantané, de taille (NxP), où Θ={θι,...,θp} 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(Θ) se décompose sous la forme : a(θp) =G(pp)Φp , avec Pp et φp représentant respectivement les paramètres non linéaires et quasi-linéaires associés à la />-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 (NxI), du bruit additif, 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 CV'7^" 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 (Λf'xNO, sera notée C^x .
FEUILLE RECTIFIEE (REGLE 91) ISA/EP
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 C2 1^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 C2^x et une étape de construction d'une fonction de coût, dite pseudo- spectre d'ordre Iq 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* (/7)" â;,vG< (/>)} de coût s'écrit sous la forme JΛ (p) = — =-7 - 7-^- , où : det{G»HG»} i\q'ψ v = ( E2^0 x) E2^x est un opérateur matriciel, dit projecteur bruit d'ordre 2q, avec E ?.2«ζ'v la matrice des vecteurs propres orthonormal i ses associés aux valeurs propres nulles de ladite matrice C2^x ;
-
Figure imgf000046_0001
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 φ.
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 œuvre une étape de minimisation de ladite fonction de coût, réalisée à partir d'une déflation à l'ordre 2q (q>\) représentative d'une récursivité d'estimation desdits P vecteurs
FEUILLE RECTIFIEE (REGLE 91) ISA/EP contenant 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 (λ≤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) ;
- calcul d'un vecteur $"" (/^) en prenant le vecteur propre correspondant à la valeur propre minimale de la matrice \ •
Figure imgf000047_0001
'
- extraction d'un vecteur <|(p) représentatif d'une estimée du vecteur de paramètres de nuisance ^(p) , à partir du vecteur
- 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 û(%p))
Figure imgf000047_0002
' Pu*s ") calcul d'une matrice
Figure imgf000047_0003
tenant compte d'un remplacement dudit vecteur a(θξ{ Λ par ledit vecteur a[θξ(pλ , et iii) réaffectation des variables suivant les fonctions suivantes :
• P ;= P-Ï ;
•
Figure imgf000047_0004
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 Â, . comprend les sous-étapes suivantes :
- extraction de M = A/9"2 vecteurs
Figure imgf000047_0005
(m) de taille (N2Xl ) ;
- transformation desdits vecteurs en M matrices BteT (m) de taille (NxN) ;
FEUILLE RECTIFIEE (REGLE 91) ISA/EP - calcul d'un vecteur propre commun aux M matrices de l'ensemble Δ^' et associé à la valeur propre la plus grande.
12. 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 œuvre une étape de construction d'un filtre FAS (pour « filtrage Alterné Séquentiel ») défini par la formule par W = [C2°*1 M®) 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.
13. Dispositif d'identification de paramètres multidimensionnels de sources d'intérêt en ce qu'il comporte un processeur adapté à mettre en œuvre 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≥\ sources présentes au sein d'un environnement multidimensionnel prédéterminé, au moyen d'une pluralité d'observations en nombre N≥\ fini, selon l'une au moins des revendications 1 à 12.
14. 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.
15. 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 œuvre 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≥\ sources d'intérêt présentes au sein d'un environnement multidimensionnel prédéterminé, au moyen d'une
FEUILLE RECTIFIEE (REGLE 91) ISA/EP pluralité d'observations en nombre N≥\ fini, selon Tune 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≥\ 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.
FEUILLE RECTIFIEE (REGLE 91) ISA/EP
PCT/EP2006/068636 2005-11-17 2006-11-17 Procede et dispositif d'identification de parametres multidimensionnels : application a la localisation et la reconstruction d'activites electriques de profondeur au moyen d'observations de surface WO2007057453A1 (fr)

Priority Applications (1)

Application Number Priority Date Filing Date Title
US12/094,078 US20090093964A1 (en) 2005-11-17 2006-11-17 Multi-dimensional parameter identification method and device: application to the location and reconstruction of deep electrical activities by means of surface observations

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
FR0511668 2005-11-17
FR0511668A FR2893434B1 (fr) 2005-11-17 2005-11-17 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

Publications (1)

Publication Number Publication Date
WO2007057453A1 true WO2007057453A1 (fr) 2007-05-24

Family

ID=36685626

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/EP2006/068636 WO2007057453A1 (fr) 2005-11-17 2006-11-17 Procede et dispositif d'identification de parametres multidimensionnels : application a la localisation et la reconstruction d'activites electriques de profondeur au moyen d'observations de surface

Country Status (3)

Country Link
US (1) US20090093964A1 (fr)
FR (1) FR2893434B1 (fr)
WO (1) WO2007057453A1 (fr)

Families Citing this family (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8321134B2 (en) * 2008-10-31 2012-11-27 Saudi Arabia Oil Company Seismic image filtering machine to generate a filtered seismic image, program products, and related methods
US20110148928A1 (en) * 2009-12-17 2011-06-23 General Electric Company System and method to correct motion in gated-pet images using non-rigid registration
WO2011103553A2 (fr) 2010-02-22 2011-08-25 Saudi Arabian Oil Company Système, machine et support de stockage lisible par ordinateur pour formation de tracé sismique amélioré à l'aide de réseau sismique virtuel
WO2021071871A1 (fr) * 2019-10-09 2021-04-15 Trustees Of Boston University Système d'électrographie utilisant des électrodes en couches en vue d'obtenir une résolution spatiale améliorée
CN111012335B (zh) * 2019-11-28 2020-10-20 燕山大学 基于非负cp分解模型的脑电意图解码方法
CN110974219A (zh) * 2019-12-20 2020-04-10 北京脑陆科技有限公司 一种基于侵入式bci的人脑意念识别系统
CN111314261B (zh) * 2020-02-24 2022-08-16 中国人民解放军国防科技大学 一种集中插入式帧同步快速盲识别方法
JP2022169141A (ja) * 2021-04-27 2022-11-09 株式会社アドバンテスト 信号ベクトル導出装置、方法、プログラム、記録媒体
CN113391362B (zh) * 2021-08-13 2021-10-29 成都理工大学 基于廊带数据约束的大地电磁剖面三维结构化反演方法
CN114459644B (zh) * 2021-12-30 2023-03-24 南京航空航天大学 基于光纤应变响应与高斯过程的起落架落震载荷辨识方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
FR2853480B1 (fr) * 2003-04-01 2005-06-17 Procede et dispositif d'identification autodidacte d'un melange sous-determine de sources au quatrieme ordre

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
ALBERA L ET AL: "Blind Identification of Overcomplete MixturEs of sources (BIOME)", LINEAR ALGEBRA AND ITS APPLICATIONS, ELSEVIER SCIENCE PUBLISHING CO., NEW YORK, NY, US, vol. 391, 1 November 2004 (2004-11-01), pages 3 - 30, XP004580007, ISSN: 0024-3795 *
BIN HE: "Modeling and Imaging of Bioelectrical Activity. Principles and Applications", 2004, KLUWER ACADEMIC/PLENUM PUBLISHERS, NEW YORK, XP002392531 *
JOHN C MOSHER ET AL: "Source Localization Using Recursively Applied and Projected (RAP) MUSIC", IEEE TRANSACTIONS ON SIGNAL PROCESSING, IEEE SERVICE CENTER, NEW YORK, NY, US, vol. 47, no. 2, February 1999 (1999-02-01), XP011058448, ISSN: 1053-587X *
SATOSHI NIIJIMA, SHOOGO UENO: "Universal Fourth Order Music Method: Incorporation of ICA into Meg Inverse Solution", PROCEEDINGS OF THE 3RD INTERNATIONAL CONFERENCE ON INDEPENDENT COMPONENT ANALYSIS AND BLIND SIGNAL SEPARATION, 9 December 2001 (2001-12-09), San Diego, California, USA, pages 639 - 644, XP002392527 *
YAO D.: "Electroencephalography inverse problem by subspace decomposition of the fourth-order cumulant matrix", SHENG WU YI XUE GONG CHENG XUE ZA ZHI / JOURNAL OF BIOMEDICAL ENGINEERING, vol. 17, no. 2, June 2000 (2000-06-01), pages 174-8, XP002392528, Retrieved from the Internet <URL:www.pubmed.gov> [retrieved on 20060727] *

Also Published As

Publication number Publication date
US20090093964A1 (en) 2009-04-09
FR2893434B1 (fr) 2008-05-09
FR2893434A1 (fr) 2007-05-18

Similar Documents

Publication Publication Date Title
WO2007057453A1 (fr) Procede et dispositif d&#39;identification de parametres multidimensionnels : application a la localisation et la reconstruction d&#39;activites electriques de profondeur au moyen d&#39;observations de surface
Albera et al. ICA-based EEG denoising: a comparative analysis of fifteen methods
Lee et al. A data-driven sparse GLM for fMRI analysis using sparse dictionary learning with MDL criterion
FR3024569A1 (fr) Procede de localisation d&#39;une activite cerebrale associee a une tache
EP2407912A2 (fr) Procédé et système de classification de signaux neuronaux, et procédé de sélection d&#39;électrodes pour commande neuronale directe
Becker et al. Multi-way space–time–wave-vector analysis for EEG source separation
Siuly et al. Identification of motor imagery tasks through CC–LR algorithm in brain computer interface
EP1906822A2 (fr) Procede et dispositif de representation d&#39;une image fonctionnelle dynamique du cerveau, par localisation et discrimination des generateurs neuroelectrioues intracerebraux et leurs applications
FR3025037A1 (fr) Procede de localisation d&#39;une activite cerebrale associee a une tache
Müller et al. Blind source separation techniques for decomposing event-related brain signals
Akojwar et al. Feature extraction of EEG signals using wavelet and principal component analysis
Becker et al. A performance study of various brain source imaging approaches
Das et al. Neuro-current response functions: A unified approach to MEG source analysis under the continuous stimuli paradigm
Wei et al. EEG-based evaluation system for motion sickness estimation
Avarvand et al. Localizing true brain interactions from EEG and MEG data with subspace methods and modified beamformers
Murzin et al. Anatomically constrained minimum variance beamforming applied to EEG
Giri et al. EEG dipole source localization in hemispherical harmonics domain
Hitziger Modeling the variability of electrical activity in the brain
Ramírez Source localization
Albera et al. Localization of spatially distributed brain sources after a tensor-based preprocessing of interictal epileptic EEG data
Becker Denoising, separation and localization of EEG sources in the context of epilepsy
EP2285273A1 (fr) Procede et systeme non invasif de detection et d&#39;evaluation de l&#39;activite electrophysiologique neuronale
Solis et al. Adaptive EEG artifact suppression using Gaussian mixture modeling
Escalona-Vargas et al. Multicompare tests of the performance of different metaheuristics in EEG dipole source localization
Ponomarev et al. Second order blind identification of event related potentials sources

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application
WWE Wipo information: entry into national phase

Ref document number: 2006830041

Country of ref document: EP

NENP Non-entry into the national phase

Ref country code: DE

WWE Wipo information: entry into national phase

Ref document number: 12094078

Country of ref document: US

122 Ep: pct application non-entry in european phase

Ref document number: 06830041

Country of ref document: EP

Kind code of ref document: A1