Method, device and equipment for denoising EEG data based on hierarchical causal graph model

Abstract

The present application discloses an EEG data denoising method, device and equipment based on a hierarchical causal graph model, the method comprising: obtaining original EEG data to be denoised and a time series corresponding to the original EEG data; segmenting the original EEG data according to the time series, obtaining a plurality of segmented data and a starting point feature vector corresponding to each segmented data; constructing a causal graph model according to the starting point feature vector corresponding to each segmented data; obtaining periodic variation information corresponding to the original EEG data according to the causal graph model; dividing the causal graph model to obtain a plurality of levels; obtaining noise information corresponding to each level according to the causal graph model and preset prior knowledge; obtaining data consistency information of a plurality of levels, and obtaining the degree of distortion of the level according to the periodic variation information and the data consistency information; constructing a target filtering algorithm according to the noise information and the degree of distortion, and filtering out the interference noise corresponding to the original EEG data at each level.

Description

Technical Field

[0001] The present application relates to the field of brain-computer interface technology, and in particular to a method, device and equipment for denoising EEG data based on a hierarchical causal graph model. Background Art

[0002] As an important means of directly measuring brain neural activity, electroencephalography (EEG) technology has broad application prospects in medical diagnosis, brain-computer interface (BCI), cognitive function analysis and other fields. Compared with other imaging techniques (such as functional magnetic resonance imaging (fMRI), EEG has the advantages of high temporal resolution, no radiation, and simple operation, and can reflect the bioelectric activity of brain neuron groups in real time. However, since the amplitude of the EEG signal itself is extremely weak (usually in the microvolt to millivolt range), it is easily interfered by various physiological and external noise sources. Therefore, removing the mixed noise components in EEG data has always been a major challenge in this field.

[0003] Traditional EEG denoising methods mainly include data-driven statistical modeling methods based on wavelet transform, independent component analysis (ICA), empirical mode decomposition (EMD), and probabilistic framework methods such as Bayesian estimation and particle filtering. These technologies perform well in simple laboratory environments, but in real scenarios, due to the lack of modeling constraints on complex noise sources and intrinsic brain mechanisms, their denoising performance is often unsatisfactory. In recent years, with the rapid development of deep learning technology, researchers have begun to try to automatically extract features from data and predict pure EEG signals based on methods such as convolutional neural networks, recursive neural networks, and generative adversarial networks, and have achieved certain results. However, these methods either lack sufficient biological explanation capabilities or cannot fully utilize prior domain knowledge, and the denoising effect has certain limitations.

[0004] Therefore, a method is urgently needed to solve at least one of the above problems. Summary of the invention

[0005] The embodiments of the present application provide a method, device and equipment for denoising EEG data based on a hierarchical causal graph model, aiming to solve the problem that existing methods either lack sufficient biological explanation capabilities or cannot fully utilize prior domain knowledge, and the denoising effect has certain limitations.

[0006] In a first aspect, an embodiment of the present application provides an EEG data denoising method based on a hierarchical causal graph model, comprising:

[0007] Acquire the original EEG data to be denoised and the time series corresponding to the original EEG data;

[0008] Segmenting the original EEG data according to the time series to obtain a plurality of segmented data and a starting point feature vector corresponding to each segmented data;

[0009] Based on the causal discovery algorithm, a causal graph model is constructed according to the starting point feature vector corresponding to each of the segmented data;

[0010] Acquiring periodic variation information corresponding to the original EEG data according to the causal graph model;

[0011] Dividing the causal graph model to obtain multiple levels; each level is used to reflect a different causal mechanism;

[0012] According to the causal graph model and preset prior knowledge, obtaining noise information corresponding to each of the levels;

[0013] Acquire data consistency information of the plurality of levels, and acquire the degree of distortion of each level according to the periodic variation information and the data consistency information;

[0014] Based on a preset causal reasoning framework, a target filtering algorithm is constructed according to the noise information and the degree of distortion, and interference noise corresponding to the original EEG data at each level is filtered out according to the target filtering algorithm.

[0015] In a second aspect, the present application also provides an EEG data denoising device, comprising:

[0016] A data acquisition module, used to acquire the original EEG data to be denoised and the time series corresponding to the original EEG data;

[0017] A data segmentation module, used to segment the original EEG data according to the time series, and obtain a plurality of segmented data and a starting point feature vector corresponding to each segmented data;

[0018] A model building module, used to build a causal graph model based on a causal discovery algorithm and according to the starting point feature vector corresponding to each segmented data;

[0019] A change acquisition module, used for acquiring periodic change information corresponding to the original EEG data according to the causal graph model;

[0020] A hierarchical division module, used to divide the causal graph model into multiple hierarchies; each hierarchy is used to reflect a different causal mechanism;

[0021] A noise acquisition module, used for acquiring noise information corresponding to each of the levels according to the causal graph model and preset prior knowledge;

[0022] A consistency acquisition module, used to acquire data consistency information of the multiple levels, and acquire the degree of distortion of each level according to the periodic variation information and the data consistency information;

[0023] The denoising completion module is used to construct a target filtering algorithm based on a preset causal reasoning framework according to the noise information and the degree of distortion, and filter out the interference noise corresponding to the original EEG data at each level according to the target filtering algorithm.

[0024] In a third aspect, the present application also provides a computer device, comprising a processor and a memory, wherein the memory is used to store a computer program, and when the computer program is executed by the processor, the method for denoising EEG data based on the hierarchical causal graph model as described in the first aspect is implemented.

[0025] In a fourth aspect, the present application also provides a computer-readable storage medium storing a computer program, which, when executed by a processor, implements the EEG data denoising method based on the hierarchical causal graph model as described in the first aspect.

[0026] Compared with the prior art, this application has at least the following beneficial effects:

[0027] 1. Achieved comprehensive and efficient multi-source noise separation and suppression

[0028] This method combines a variety of theories and technologies to implement comprehensive and multi-level identification and removal of various noise sources. Through hierarchical causal modeling and prior knowledge, this method can estimate and quantify the different types of noise and their impact on each level, laying the foundation for noise separation. In layered and block filtering and local network asynchronous noise separation, targeted and differentiated processing of noise in different levels and different network modules is achieved, effectively overcoming the defects of the traditional "one-size-fits-all" filtering method.

[0029] This method can accurately separate and remove various common physiological noises (such as electrooculography and electromyography), external noises (such as electromagnetic interference) and internal noise sources (such as distraction), ensuring that the output EEG data x^*(t) is purified to the maximum extent. At the same time, through innovative leading technologies such as hierarchical distortion measurement and network synchronization analysis, this method can automatically identify the scope of noise influence and avoid unnecessary damage to uncontaminated valid signals. This focused and intelligent hierarchical and block denoising paradigm can maximize the elimination of interference noise while retaining the integrity and signal details of the original EEG data to the greatest extent, achieving a win-win situation of "denoising and enhancement".

[0030] 2. Improved signal robustness and application reliability

[0031] Since this method can effectively suppress the superposition of multi-source noise, it will inevitably greatly improve the signal quality and robustness of the output data x^*(t). This method comprehensively evaluates the degree of synchronization between different brain regions through multi-dimensional comprehensive analysis of time-frequency amplitude, phase, etc., quickly identifies areas that have lost synchronization due to noise interference, and focuses on repairing these "abnormal" areas. This synchronization-based analysis paradigm can ensure the efficient transmission of neural oscillation signals between levels and networks, and maintain the normal coordination between different functional modules of the brain.

[0032] In addition, through the "no residue processing" finishing process (wavelet filtering, sparse decomposition, etc.) of multi-technology integration, this method can also implement a "sweep-up" type of deep denoising for residual wavelet noise and transient transient noise, further improving the overall quality and robustness of the pure signal x^(t). High-quality x^(t) provides a reliable basic data guarantee for subsequent analysis and application, and can significantly improve the credibility and accuracy of EEG data analysis and interpretation.

[0033] It should be understood that the foregoing general description and the following detailed description are exemplary and explanatory only and are not restrictive of the present application. BRIEF DESCRIPTION OF THE DRAWINGS

[0034] Figure 1 A schematic diagram of a process of an EEG data denoising method based on a hierarchical causal graph model according to an embodiment of the present application;

[0035] Figure 2 This is a schematic diagram of the structure of an EEG data denoising device shown in an embodiment of the present application;

[0036] Figure 3 A schematic diagram of the structure of a computer device shown in an embodiment of the present application. DETAILED DESCRIPTION

[0037] In the following description, specific details such as specific system structures, technologies, etc. are provided for the purpose of illustration rather than limitation, so as to provide a thorough understanding of the embodiments of the present application. However, it should be clear to those skilled in the art that the present application may also be implemented in other embodiments without these specific details. In other cases, detailed descriptions of well-known systems, devices, circuits, and methods are omitted to prevent unnecessary details from obstructing the description of the present application.

[0038] It should be understood that when used in the present specification and the appended claims, the term "comprising" indicates the presence of described features, wholes, steps, operations, elements and / or components, but does not exclude the presence or addition of one or more other features, wholes, steps, operations, elements, components and / or combinations thereof.

[0039] It should also be understood that the term “and / or” used in the specification and appended claims refers to any and all possible combinations of one or more of the associated listed items, and includes these combinations.

[0040] As used in the specification and appended claims of this application, the term "if" can be interpreted as "when" or "uponce" or "in response to determining" or "in response to detecting", depending on the context. Similarly, the phrase "if it is determined" or "if [described condition or event] is detected" can be interpreted as meaning "uponce it is determined" or "in response to determining" or "uponce [described condition or event] is detected" or "in response to detecting [described condition or event]", depending on the context.

[0041] In addition, in the description of the present application specification and the appended claims, the terms "first", "second", "third", etc. are only used to distinguish the descriptions and cannot be understood as indicating or implying relative importance.

[0042] References to "one embodiment" or "some embodiments" etc. described in the specification of this application mean that one or more embodiments of the present application include specific features, structures or characteristics described in conjunction with the embodiment. Therefore, the statements "in one embodiment", "in some embodiments", "in some other embodiments", "in some other embodiments", etc. that appear in different places in this specification do not necessarily refer to the same embodiment, but mean "one or more but not all embodiments", unless otherwise specifically emphasized in other ways. The terms "including", "comprising", "having" and their variations all mean "including but not limited to", unless otherwise specifically emphasized in other ways.

[0043] The technical solution of the embodiment of the present application is introduced below.

[0044] As an important means of directly measuring brain neural activity, electroencephalography (EEG) technology has broad application prospects in medical diagnosis, brain-computer interface (BCI), cognitive function analysis and other fields. Compared with other imaging techniques (such as functional magnetic resonance imaging (fMRI), EEG has the advantages of high temporal resolution, no radiation, and simple operation, and can reflect the bioelectric activity of brain neuron groups in real time. However, since the amplitude of the EEG signal itself is extremely weak (usually in the microvolt to millivolt range), it is easily interfered by various physiological and external noise sources. Therefore, removing the mixed noise components in EEG data has always been a major challenge in this field.

[0045] Traditional EEG denoising methods mainly include data-driven statistical modeling methods based on wavelet transform, independent component analysis (ICA), empirical mode decomposition (EMD), and probabilistic framework methods such as Bayesian estimation and particle filtering. These technologies perform well in simple laboratory environments, but in real scenarios, due to the lack of modeling constraints on complex noise sources and intrinsic brain mechanisms, their denoising performance is often unsatisfactory. In recent years, with the rapid development of deep learning technology, researchers have begun to try to automatically extract features from data and predict pure EEG signals based on methods such as convolutional neural networks, recursive neural networks, and generative adversarial networks, and have achieved certain results. However, these methods either lack sufficient biological explanation capabilities or cannot fully utilize prior domain knowledge, and the denoising effect has certain limitations.

[0046] Please refer to Figure 1 , Figure 1 A flowchart of a method for denoising EEG data based on a hierarchical causal graph model provided in an embodiment of the present application. The method for denoising EEG data based on a hierarchical causal graph model in an embodiment of the present application can be applied to computer devices, including but not limited to smartphones, laptops, tablet computers, desktop computers, physical servers, cloud servers and other devices. Figure 1 As shown, the EEG data denoising method based on the hierarchical causal graph model of this embodiment includes steps S101 to S108, which are described in detail as follows:

[0047] Step S101, obtaining original EEG data to be denoised and a time series corresponding to the original EEG data.

[0048] Specifically, a portable EEG headband device is used to collect EEG signals. The headband has multiple built-in sensing electrodes, usually dry electrodes or wet electrodes. Dry electrodes do not need to use conductive paste, and can detect EEG signals as long as they are in good contact with the scalp, which is convenient for ordinary users to operate by themselves. Wet electrodes need to use conductive paste as a medium. Although the operation is a bit cumbersome, the signal quality is better. During the collection process, the user first wears the EEG headband, and adjusts the position and tightness of the headband according to the prompts so that each sensing electrode can be close to the scalp.

[0049] Next, the user needs to connect and start the EEG acquisition application through a mobile device (such as a mobile phone or tablet). The application will automatically detect the connection status of the headband and the electrode signal strength, and prompt the user to make necessary adjustments based on the feedback, such as re-adjusting the headband position, applying more conductive paste, etc. At the same time, the application will also guide the user to relax physically and mentally, and try to avoid muscle movement and distraction to reduce interference from electromyography and other physiological noise. When all indicators meet the requirements, the application will prompt the user to enter the formal EEG acquisition stage.

[0050] During the EEG acquisition phase, users may need to perform specific cognitive or motor tasks depending on the needs of the specific application scenario. For example, in a motor imagery application, users need to imagine boxing or kicking a football according to the prompts, and the EEG device will synchronously record the changes in the brain's electrical activity at this time. In an attention detection application, users need to focus on the stimuli on the screen, and the device will capture changes in attention levels. In addition, there are some application scenarios where users do not need to perform specific tasks, but only need to relax quietly, and the device will continue to record baseline EEG activity at this time.

[0051] Regardless of the application scenario, the EEG acquisition device will record the weak potential changes detected by the sensor electrodes in real time. After amplification, filtering and digital processing, these analog electrical signals will be stored in the device memory or cloud in the order of time series and channel number to form a multi-channel, large-capacity raw data file. The European Data Format (EDF) or the biological semifEDF format is usually used to store EEG data, including metadata header files and binary data files. Taking mobile terminal acquisition as an example, this application can set the sampling rate to 256Hz, 16-bit quantization accuracy, 0.1-100Hz bandpass filtering, collect 64-channel EEG data for 10 minutes, and generate a raw data file of about 75MB.

[0052] Step S102, segmenting the original EEG data according to the time series, obtaining a plurality of segmented data and a starting point feature vector corresponding to each segmented data.

[0053] Specifically, the original EEG data is recorded continuously in time series. In order to facilitate subsequent causal modeling and analysis, this application needs to divide it into multiple time segments. The basic unit of segmentation is the EEG epoch, which refers to a continuous data segmented from the original time series, usually 1-2 seconds in length. Dividing the original data into multiple epochs can effectively capture the instantaneous changes and local features in the EEG signal.

[0054] Taking the original EEG data stored in EDF format as an example, first read the sampling rate fs, number of channels N, data length L and other parameters from the header file (.HDR). Assume that the sampling rate fs = 512Hz, the number of channels N = 64, and the data length L = 307200 sample points (corresponding to 10 minutes). This application divides the time series into M = 600 epochs, that is, each epoch is 512 sample points (corresponding to 1 second). For the i-th epoch (1≤i≤M), its start time is (i-1)512 / fs, and its end time is i512 / fs. In the process of equal division, if the subject has large movements or other abnormalities, the time period can be eliminated according to the marking information to avoid introducing too much noise.

[0055] Next, for each time segment epoch, this application needs to extract its starting point feature vector. This feature vector reflects the local characteristics of the EEG signal such as amplitude and phase at the beginning of the epoch, which will provide important clues for subsequent causal discovery. The specific method is: take the first data point of the epoch as the starting point, take the starting point as the origin, take the subsequent n data points (n is a preset constant, such as n=10), and form a (n+1)-dimensional vector of these n+1 data points in time series, which is the starting point feature vector of the epoch.

[0056] Suppose the j-th channel data of the i-th epoch is xi,j=[xi,j(0),xi,j(1),...,xi,j(511)]^T, where xi,j(k) is the sampling value at the k-th time point. Then the starting point feature vector of the channel is:

[0057] vi,j=[xi,j(0),xi,j(1),...,xi,j(n)]^T.

[0058] Concatenate the starting point feature vectors of all N channels to form a total feature vector of (n+1)*N dimensions:

[0059] Vi=[vi,1^T,vi,2^T,...,vi,N^T]^T.

[0060] This total feature vector Vi is the starting point feature vector of the i-th epoch, with a dimension of (n+1)*N. Taking the previous example, assuming n=10, Vi is a 74-dimensional vector. In this way, the application converts the original time series EEG data into the form of a local feature vector.

[0061] Step S103, based on the causal discovery algorithm, construct a causal graph model according to the starting point feature vector corresponding to each of the segmented data.

[0062] Specifically, this application proposes an innovative causal discovery algorithm - the FCI-GNN algorithm based on graph neural network and Monte Carlo sampling, which is used to learn the causal graph model from the starting point feature vector of EEG data. The core idea is: first, use the constraint-based FCI algorithm to recover an approximate skeleton graph from the observed data, and then use the skeleton graph as input to repair and orient the edges of the skeleton graph through an end-to-end graph neural network model, and finally output a complete directed acyclic graph (DAG) as the causal graph model.

[0063] In some embodiments, the causal discovery algorithm is based on constructing a causal graph model according to the starting point feature vector corresponding to each of the segmented data, including: based on the skeleton recovery algorithm, obtaining an undirected skeleton graph in each of the starting point feature vectors; inputting the undirected skeleton graph into a graph neural network encoder, and outputting a representation vector corresponding to the undirected skeleton graph; sampling each of the representation vectors according to a Monte Carlo graph sampling strategy to obtain multiple candidate directed acyclic graph samples; obtaining a candidate graph score corresponding to each candidate directed acyclic graph sample, and determining a target sample from the multiple candidate directed acyclic graph samples according to the candidate graph score; inputting the target sample into a preset graph neural network decoder, and outputting the causal graph model.

[0064] Specifically, the algorithm includes the following key steps:

[0065] (1) Skeleton recovery. This application first runs the classic FCI skeleton recovery algorithm to obtain an undirected skeleton graph Gi from the starting point feature vector Vi (dimension is (n+1)*N) of each epoch. The skeleton graph reflects the disordered correlation between the observed variables.

[0066] Exemplarily, the obtaining of the undirected skeleton graph in each of the starting point feature vectors includes: forming a completely undirected graph based on a plurality of the starting point feature vectors; performing a conditional independence test on the completely undirected graph, deleting edges in the completely undirected graph that do not satisfy a preset separation criterion, and obtaining the undirected skeleton graph.

[0067] Specifically, the skeleton recovery steps of the FCI algorithm include: first, constructing the (n+1)*N features in Vi into a completely undirected graph, with each feature as a node and an undirected edge between any two nodes; then, based on conditional independence tests, such as the chi-square test and mutual information test, gradually delete those edges that do not meet the d-separation criterion until the largest equivalence class that meets the conditional independence is found.

[0068] Taking the i-th epoch as an example, let Vi=[vi,1^T, vi,2^T, ..., vi,N^T]^T, where vi,j (j=1,...,N) is the starting point feature vector of the j-th channel, then the initial completely undirected graph contains (n+1)*N nodes. After conditional independence test and edge deletion, this application obtains an undirected skeleton graph Gi, in which any two connected nodes indicate that there is a correlation between the corresponding two features, but the nature of the correlation (causality or common cause) is not clear.

[0069] (2) Graph encoding. The skeleton graph Gi is input into the encoder of the graph neural network, and the graph data is encoded using models such as the graph convolutional neural network (GCN) to learn the representation vector of each node. In this process, GCN can well capture the graph structure information and node attributes, providing sufficient prior knowledge for the next sampling process.

[0070] Specifically, the present application represents the undirected skeleton graph Gi as G=(V,E), where V is a node set and E is an undirected edge set. Each node v∈V corresponds to a node embedding vector h_v, which is used to encode the attribute information of the node. The core idea of ​​GCN is to update the node representation through information transfer between nodes, making full use of the prior knowledge of the graph structure. In each layer of GCN, the embedding vector of node v is obtained by combining the embedding vectors of its neighboring nodes and its own embedding vector:

[0071] h_v^(l+1) = σ(W^(l)·COMBINE({h_u^(l),∀u∈N(v)})).

[0072] Among them, N(v) is the set of neighbor nodes of node v, COMBINE is a symmetric function (such as summation, averaging, etc.) used to summarize the embedding vectors of neighbor nodes; W^(l) is the trainable weight matrix of this layer; σ is a nonlinear activation function. Through the multi-layer GCN structure, this application finally obtains the final embedding vector vector h_v^(L) of each node v as the representation vector of the node.

[0073] In this encoding process, GCN can well capture graph structure information (i.e., the topological relationship between nodes) and node attributes (such as the numerical value of the corresponding feature), providing sufficient prior knowledge for the subsequent sampling process.

[0074] (3) Monte Carlo graph sampling. Based on the node representation learned in step (2), this application uses a Monte Carlo sampling strategy to perform random walks in the candidate DAG space to generate a large number of candidate directed acyclic graph samples. This sampling process makes full use of the prior knowledge extracted by the graph neural network encoder and can explore possible DAG spaces more efficiently.

[0075] Specifically, this application uses the node representation h_v^(L) obtained by GCN encoding as input, and designs a scoring function f(G|h_v^(L)) to measure the probability (or likelihood) score of the directed acyclic graph G under a given node representation. Then, this application uses a sampling strategy based on Markov Chain Monte Carlo (MCMC) to perform random walks in all possible DAG spaces and generate a large number of candidate DAG samples.

[0076] Taking the Metropolis-Hastings algorithm as an example, in each sampling step, this application starts from the current DAG sample G_t and obtains a new candidate DAG G_new through certain graph editing operations (such as adding / deleting edges, reversing edge directions, etc.). Then calculate the acceptance rate of G_new relative to G_t:

[0077] α = min(1, f(G_new|h_v^(L)) / f(G_t|h_v^(L))).

[0078] Next, this application accepts G_new as the next sample G_t+1 with probability α, or rejects G_new with probability 1-α and retains G_t. Through multiple iterations, a large number of high-quality DAG samples can be obtained.

[0079] This sampling process makes full use of the prior knowledge extracted by the GCN encoder, allowing this application to explore the entire DAG space more efficiently without falling into a local optimal solution. At the same time, due to the theoretical guarantee of MCMC, the DAG sample set finally obtained by this application obeys the target distribution f(G|h_v^(L)).

[0080] (4) Candidate graph scoring. For each sampled candidate DAG, this application decodes it into a standard causal graph representation and scores it using a scoring function (such as BDeu score) to measure its consistency with the observed data.

[0081] Specifically, for the i-th epoch, this application obtains a series of candidate DAG samples {G_1, G_2, ..., G_M}. It should be noted that these DAG samples do not directly reflect the causal structure of the observed data Vi, because during the encoding and sampling process, this application only uses the topological structure information of the skeleton graph Gi, without considering the observed data Vi itself.

[0082] Therefore, this application needs to decode each DAG G_m into the corresponding standard causal graph representation CGM(G_m), and perform parameter learning on Vi according to CGM(G_m) to obtain the likelihood score LL(Vi|CGM(G_m)) given Vi. This application uses this likelihood score as the score of DAG G_m to measure its fit with the observed data. Commonly used parameter learning and scoring methods include BDeu score, BIC score, etc.

[0083] (5) Graph decoding and output. The candidate DAG with the highest score is taken as the final output, which is the causal graph model G of this epoch. Through the decoder of the graph neural network, this application can decode the DAG into a standard causal graph representation, intuitively presenting the causal relationship between variables.

[0084] Specifically, this application selects the DAG with the highest score from the candidate DAG sample set {G_1, G_2, ..., G_M}, denoted as G*, that is:

[0085] G* = arg max_{G_m} LL(Vi|CGM(G_m)).

[0086] Decode G into a causal graph representation CGM(G), and this application obtains the final causal graph model of the i-th epoch:

[0087] G_i = CGM(G*).

[0088] The causal graph model Gi is a directed acyclic graph, where each node represents a feature dimension in Vi, and the directed edge X→Y indicates that X is the direct causal cause of Y.

[0089] The decoding process is essentially to convert the abstract DAG structure G into an explicit causal graph representation CGM(G). This application can design a structured decoder that takes each node in G and its neighbor node information as input and outputs the causal relationship representation of the node in CGM(G).

[0090] For example, a sequence-to-sequence (Seq2Seq) model can be used, where the encoder encodes each node in G* and its neighbor topology, and the decoder generates a causal statement for the node in CGM, such as "X→Y", "X←Y", "X||Y" (indicating that X and Y have no causal relationship), etc. Through the attention mechanism, the decoder can focus on different neighbors of the input node and flexibly generate causal statements for different topological structures.

[0091] The above-mentioned FCI-GNN algorithm makes full use of the constraint advantages of the classic FCI algorithm and the advantages of graph neural network modeling nonlinearity and capturing hidden variables, and can more accurately learn the inherent causal mechanism from complex EEG data.

[0092] Taking a simplified 10-dimensional example to illustrate, assume that the FCI-GNN algorithm finally outputs the following causal graph model:

[0093] G = ({X1, X2, X3, ..., X10}, {X1->X3, X2->X3, X4->X7, X5->X7, X7->X8,

[0094] This directed acyclic graph G contains 10 nodes, each node Xi (1≤i≤10) represents the i-th dimension feature, that is, the starting point feature of the i-th EEG channel in the epoch. There are 9 directed edges in the graph G, and each directed edge Xi->Xj means that Xi is the direct cause of Xj, that is, Xi has a direct causal effect on Xj.

[0095] Among them, the edges X1->X3 and X2->X3 indicate that X1 and X2 are both direct causes of X3, which may reflect that the two brain regions corresponding to X1 and X2 jointly affect the electrical activity of the X3 brain region; the edges X4->X7 and X5->X7 indicate that X4 and X5 are the causes of X7, that is, the activities of the X4 and X5 brain regions lead to changes in the activity of the X7 brain region; X7->X8 indicates that X7 is the cause of X8, reflecting the existence of a causal transmission relationship between the X7 and X8 brain regions; X6->X9 and X9->X10 reveal the hierarchical causal links between the three brain regions X6, X9, and X10.

[0096] It is worth noting that there is a special edge X5->X1 in the causal graph, indicating that X5 has a direct causal effect on X1. However, according to prior knowledge, there should not be a direct causal effect between the brain regions corresponding to X1 and X5, so this edge X5->X1 is likely to be caused by the influence of hidden variables. In other words, there is an unobserved hidden brain region that affects the electrical activities of both X1 and X5, resulting in a direct causal relationship between X1 and X5 in the observed data. The FCI-GNN algorithm can effectively capture this hidden variable influence through the graph neural network and reflect it in the final causal graph model.

[0097] Step S104, obtaining periodic variation information corresponding to the original EEG data according to the causal graph model.

[0098] Specifically, this application uses an innovative graph convolutional neural network-based periodicity detection algorithm (PGCN) to identify and quantify the periodic components of each epoch based on the causal graph model constructed in steps. Under the guidance of the causal graph model, the PGCN algorithm can adaptively capture the inherent periodic oscillations of EEG signals and quantify the degree of variation of each oscillation cycle, providing important support for subsequent hierarchical modeling and noise separation.

[0099] In some embodiments, obtaining the periodic variation information corresponding to the original EEG data according to the causal graph model includes: identifying the periodic embedding vector corresponding to each of the causal graph models according to a periodic detection algorithm based on a graph convolutional neural network; decoding the periodic embedding vector to obtain periodic parameters; and generating the periodic variation information according to the periodic parameters.

[0100] Specifically, for the causal graph model G_i of the i-th epoch, this application inputs it into the PGCN algorithm for processing. The core of the PGCN algorithm is to design a new type of graph convolution layer that can learn periodic representations on the graph structure. This application defines G_i=({X_1, X_2, ..., X_N}, E) as the input causal graph, where X_j (j=1,...,N) is the j-th feature dimension, N=(n+1)*M, and E is a set of directed edges. For each node X_j, this application constructs its k-neighbor set:

[0101] N_k(X_j) = {X_t | d(X_j, X_t) ≤ k, X_t ∈ V}.

[0102] Where d(X_j, X_t) represents the shortest path length between X_j and X_t. Then the k-neighbor embedding vector of X_j can be defined as:

[0103] h_j^k = ψ(AGG({X_t, ∀X_t ∈ N_k(X_j)})).

[0104] Here AGG is an aggregation function (such as summation, mean, etc.) of a node and its k-neighbors, and ψ is a trainable projection function used to map the aggregated embedding vector to a vector space of fixed dimension. In this way, the application can capture the comprehensive information of node X_j and its k-neighbors, that is, the topological structure and attribute characteristics of the node in the area with a radius of k.

[0105] Furthermore, the present application stacks the k-neighbor embedding vectors into multi-scale embedding vectors:

[0106] H_j = [h_j^1, h_j^2, ..., h_j^K].

[0107] Where K is the maximum neighbor radius, H_j ∈ R^(K*d) is the multi-scale embedding vector of node X_j, which can simultaneously encode graph topology and node attribute information at different scales.

[0108] In PGCN, this application designs a new convolutional layer with input H_j and output as the periodic embedding vector p_j of the node: p_j = ReLU(W_p * H_j + b_p).

[0109] Where W_p and b_p are trainable parameters, and ReLU is a nonlinear activation function. This new convolutional layer can adaptively learn periodic representations from the input multi-scale embedding vectors and ultimately output a periodic embedding vector p_j.

[0110] Through the calculation of the above PGCN layer, the present application can obtain the periodic embedding vector p_j of each node X_j. Next, the present application designs a novel periodic decoder to decode p_j into the periodic parameters of the node, including:

[0111] Amplitude A_j: It represents the amplitude variation range of the EEG channel corresponding to this feature dimension.

[0112] Period T_j: Characterizes the dominant oscillation period of this characteristic dimension.

[0113] Phase φ_j: Characterizes the phase offset of this feature dimension.

[0114] Specifically, the decoder is a feed-forward neural network with input p_j and output (A_j, T_j, φ_j). Through end-to-end training, the decoder is able to learn to map from periodic embedding vectors to explicit periodic parameters.

[0115] Taking the i-th epoch as an example, after processing by PGCN and the periodic decoder, this application can obtain the periodic parameters of each feature dimension in epoch i:

[0116] {(A_1, T_1, φ_1), (A_2, T_2, φ_2), ..., (A_N, T_N, φ_N)}.

[0117] Where N=(n+1)*M, each (A_j, T_j, φ_j) represents the amplitude, period and phase of the jth feature dimension respectively. This set of periodic parameters fully describes the periodic components and their changes of the EEG signal in epoch i.

[0118] For example, for the alpha (8-13Hz) and beta (13-30Hz) frequency bands, if this application finds that the dominant period T_j of a certain feature dimension falls within these two frequency bands, it means that the feature dimension (corresponding EEG channel) has alpha or beta oscillation activity in this epoch; the amplitude A_j reflects the strength of the oscillation, and the phase φ_j describes the oscillation time offset of the channel relative to other channels.

[0119] Step S105, dividing the causal graph model to obtain multiple levels; each level is used to reflect a different causal mechanism.

[0120] Specifically, the present application performs topological sorting on the original causal graph model G=({X1, X2, ..., XN}, E) to obtain an ordered node sequence. Then, based on the shortest path distance between nodes and prior knowledge, the present application divides the nodes into L levels {L1, L2, ..., LL}. Generally speaking, the node Xj in the low-level Li is connected to more direct causes and direct effects, reflecting a more direct causal influence; while the node Xp in the high-level Lk (k>i) is connected to more indirect causes and indirect effects, reflecting a deep multi-step causal relationship.

[0121] In some embodiments, dividing the causal graph model to obtain multiple levels includes: topologically sorting the causal graph model to obtain an ordered node sequence; obtaining node path information corresponding to the ordered node sequence; and dividing the ordered node sequence into multiple levels according to the node path information.

[0122] The specific level division algorithm is as follows: Assuming that the ordered node sequence [X1, X2, ..., XN] has been obtained, let cnt[i]=0 (1≤i≤N) to indicate that the level number of the i-th node is initialized to 0, and maxlevel=0 is the maximum level number. This application traverses the ordered node sequence, and for each node Xi, according to the level of its direct cause Direct_Causes(Xi) and direct effect Direct_Effects(Xi), the level cnt[i] of the node is determined as: cnt[i] = 1 + max{cnt[j] |Xj ∈ Direct_Causes(Xi)} where the direct cause Direct_Causes(Xi) and direct effect Direct_Effects(Xi) can be directly obtained from the original causal graph G. At the same time, this application maintains the global maximum level maxlevel to ensure that cnt[i] does not exceed the current maxlevel, otherwise update maxlevel. Through this dynamic programming method, this application finally obtains the level cnt[i] of each node Xi, and the number of the bottom node must be smaller than the high-level node. According to the level number, this application classifies all nodes into L1, L2, ..., LmaxLevel, thus completing the level division.

[0123] For example, suppose the original causal graph G contains 10 nodes, in topological order [X4, X6, X8, X3, X1, X7,X5, X9, X2, X10], and has the following causal edges: Direct_Causes(X4)={}, Direct_Effects(X4)={X6, X8} Direct_Causes(X6)={X4}, Direct_Effects(X6)={X3, X8} ... Direct_Causes(X10)={X5, X9}, Direct_Effects(X10)={} Then after hierarchical division, the node sets of each layer are: L1={X4}, L2={X6, X8}, L3={X3, X1, X7}, L4={X5, X9, X2}, L5={X10}

[0124] It can be seen that in the bottom layer L1, there is only node X4, which has no direct cause and reflects the initial source of EEG activity; then nodes in L2 and L3 layers such as X6, X8, X3, X1, X7, etc. have direct or multi-step causal relationships with the nodes in L1, reflecting a higher-level regulatory mechanism; finally, nodes in L4 and L5 layers such as X2, X5, X9, X10, etc. reflect a deeper level of complex causal network.

[0125] Step S106, obtaining noise information corresponding to each of the levels according to the causal graph model and preset prior knowledge.

[0126] Specifically, based on the existing physiological and neuroscience knowledge, this application summarizes the main noise sources that may exist at different levels and their impact patterns to constitute noise information. In the lowest level (L1), the noise mainly comes from the sensor acquisition error and low-level physiological interference, such as high-frequency noise caused by poor electrode contact, myoelectric noise generated by muscle movement, etc. These bottom-level noises usually show the characteristics of Gaussian white noise or impulse noise, and the impact range is relatively local. In the middle level (such as L2 and L3), in addition to the bottom-level noise, there are also some low-frequency noises inherent in the central nervous system, such as the influence of biological currents such as respiratory filtering and blood circulation. These noises often show low-frequency drift or oscillation patterns and have a wide range of impact. In the top level (such as L4 and L5), the noise mainly comes from the inherent disorders of the brain's advanced cognitive functions, such as distraction, emotional fluctuations, etc., which manifest as nonlinear trough noise or burst noise, with the widest impact range. In addition, external electromagnetic interference will also become a common noise source for all levels.

[0127] In some embodiments, obtaining the noise information corresponding to each of the levels includes: obtaining the connection topology attributes corresponding to each of the levels; obtaining the noise propagation path and influence range according to the connection topology attributes; and generating the noise information according to the noise propagation path and influence range.

[0128] This application analyzes the connection topology properties of nodes (brain regions) in each layer based on the layered causal graph model in the step, and identifies possible noise propagation paths and impact ranges. Specifically, for the node Xi in the kth layer Lk, this application defines the following three indicators:

[0129] 1. Outdegree (Xi) is the size of Xi’s direct influence set, that is, the number of other nodes that Xi has direct causal influence on.

[0130] 2. Indegree(Xi) is the number of sources directly affected by Xi.

[0131] 3. Independent_Paths(Xi) is the number of different paths from the root node of the hierarchy to Xi.

[0132] It can be expected that a node Xi with a larger out-degree and number of independent paths is more likely to propagate noise to downstream nodes because it has multiple influencing paths and direct impacts on the downstream. A node Xj with a larger in-degree is more susceptible to the superposition of multiple noise sources upstream. Therefore, this application can regard Outdegree(Xi) as a quantitative indicator of the degree of noise impact on the downstream, and Indegree(Xi) as a quantitative indicator of the cumulative propagation of noise. Independent_Paths(Xi) measures the number of independent paths that noise may reach the node. By comparing the three indicators of nodes within the hierarchy, this application can roughly estimate the scope and intensity of the impact of different noise sources at the hierarchy.

[0133] Specifically, the present application can define the noise impact index NLk of the level Lk as follows: NLk = Σ(Xi∈Lk)[w1Outdegree(Xi) + w2Indegree(Xi) + w3*Independent_Paths(Xi)] where w1, w2, and w3 are adjustable weight coefficients used to balance the relative importance of the three indicators. The larger the NLk, the more serious the level is affected by the noise. The present application can also further compare the noise impact indexes of different levels to determine the main level of noise. For example, for the bottom layer L1, the outdegree term Outdegree(Xi) contributes the most because this level is mainly responsible for the generation of the initial input; while for the high-level L4 and L5, the indegree term Indegree(Xi) and the path term Independent_Paths(Xi) contribute more because the noise is mainly accumulated and propagated to the control core through multi-path pathways.

[0134] In addition, it is necessary to combine the prior functional properties of each level node (brain area) and the known noise pattern to make a more precise estimate of the possible noise sources and their impact. For example, through the cortical division map, this application knows that the frontal cortex is mainly responsible for cognitive control functions and is easily affected by internal noise such as distraction, while the occipital cortex is more susceptible to visual stimulation and external physical noise. For another example, based on the known event-related potential (ERP) components, this application can attribute high-level cognition-related ERP components such as P300 to high-level attention control centers, and attribute early visual-related components such as N170 to primary visual processing areas, thereby pre-identifying possible noise sources.

[0135] Finally, based on the noise impact index and prior functional attributes, this application constructs a hierarchical noise impact model to quantitatively describe the degree to which each level is affected by different noise sources. Taking the i-th epoch as an example, its hierarchical noise impact model NMi is as follows: NMi = {(Lk, NSk) | 1 ≤ k ≤ maxlevel} where NSk is the noise source distribution of the k-th layer Lk, denoted as: NSk= {(n1, p1), (n2, p2), ... , (nP, pP)} Each tuple (np, pp) indicates that the probability of noise type np affecting the k-th layer is pp, p1+p2+...+pP=1. Noise type np includes physiological noise (such as electromyography, electrooculography, etc.), external noise (such as electromagnetic interference, environmental noise, etc.) and intrinsic noise (distraction, nonlinear coupling, etc.). For the i-th epoch, this application first estimates the noise type set NS={n1, n2, ..., nP} and its prior probability P(np) according to its experimental conditions and environmental parameters, and then calculates the normalized impact probability of each noise type at each level by combining factors such as the hierarchical noise impact index NLk: pp(k) = P(np) * NLk / Σ(1≤j≤maxlevel) NLj to obtain a complete hierarchical noise impact model NMi. This model comprehensively describes the various noise sources that may exist at different levels and their impact ratios.

[0136] Step S107 , obtaining data consistency information of the multiple levels, and obtaining the distortion degree of each level according to the periodic variation information and the data consistency information.

[0137] First, this application defines the intra-level data consistency indicator IntraLevelCons(Lk) to quantify the consistency or similarity of the internal data points of the kth layer Lk. This indicator is calculated based on the differences in the periodic parameters of each channel and statistics such as mutual information. Specifically, assuming that Lk contains Mk nodes, that is, Lk={X1,X2,...,XMk}, each node Xi corresponds to the time series data xi and its periodic parameters (Ai,Ti,φi), then the internal consistency of the kth layer is:

[0138] IntraLevelCons(Lk) = Σ(1≤i,j≤Mk) ω1dA(Ai,Aj) + ω2dT(Ti,Tj) + ω3dφ(φi,φj) + ω4dMI(xi,xj).

[0139] Where dA(Ai,Aj) is the amplitude difference between nodes Xi and Xj, which can be calculated using |Ai-Aj| or other amplitude distance metrics; dT(Ti,Tj) is the period difference, which can be measured using |Ti-Tj| or other period distance metrics; dφ(φi,φj) is the phase difference, which can be calculated using min(|φi-φj|,2π-|φi-φj|); dMI(xi,xj) is the channel mutual information distance, which reflects the statistical correlation between xi and xj. ω1~ω4 are adjustable weights used to balance the importance of different indicators. The smaller the value of IntraLevelCons(Lk), the higher the data consistency within the kth layer and the lower the distortion.

[0140] Secondly, this application defines the inter-level data consistency indicator InterLevelCons(Lk, Lt) to quantify the consistency of data between Lk and Lt at different levels. This indicator focuses on the data coupling and causal transmission relationship between levels. Assuming that Lk and Lt contain Mk and Mt nodes respectively, the steps to calculate the consistency between levels are as follows:

[0141] 1. Based on the hierarchical causal model in step 1, determine the causal path set P between Lk and Lt;

[0142] 2. For each path p∈P, calculate its path length len(p);

[0143] 3. For each pair of adjacent nodes (Xi, Xj) on path p, calculate their periodic component difference d(Xi, Xj), for example, d(Xi, Xj) = ω1dA(Ai, Aj) + ω2dT(Ti, Tj) + ω3*dφ(φi, φj);

[0144] 4. Define the total difference measure of path p as D(p)=Σ(Xi,Xj∈p)d(Xi,Xj)f(len(p)), where f(len(p)) is the path length loss factor (e.g., f(len(p))=exp(-βlen(p))), which is used to penalize overly long paths;

[0145] 5. Calculate the average difference metric D̅ = (Σ(p∈P)D(p)) / |P| of all paths p, as the negative value of the inter-level consistency indicator InterLevelCons(Lk,Lt).

[0146] The larger the value of InterLevelCons(Lk, Lt), the higher the data coupling consistency between Lk and Lt at different levels, which is more conducive to the efficient transmission of noise-free signals between levels. On the contrary, when the value of this indicator is low, it means that noise has accumulated and distorted the data on the causal propagation path between levels, and effective separation is needed.

[0147] Finally, this application combines the intra-level consistency IntraLevelCons(Lk) and the inter-level consistency InterLevelCons(Lk,Lt) to obtain the overall distortion metric Distortion(Lk) of the kth layer Lk:

[0148] Distortion(Lk) = [1 - IntraLevelCons(Lk)] + Σ(t≠k) [1 -InterLevelCons(Lk,Lt)].

[0149] The first item reflects the degree of data distortion within Lk, and the second item measures the cumulative data distortion between Lk and all other layers. The larger the value of Distortion(Lk), the more serious the data distortion of the kth layer is, and strong denoising is required. This application can use Distortion(Lk) as a noise separation criterion. For layers with larger Distortion(Lk) values, stronger filtering intensity will be applied to ensure signal quality.

[0150] In some embodiments, before constructing the target filtering algorithm according to the noise information and the degree of distortion, it also includes: obtaining frequency domain feature information corresponding to the original EEG data; obtaining brain region synchronization evaluation results according to the frequency domain feature information and periodic variation information, so as to construct a target filtering algorithm according to the noise information, the degree of distortion, and the brain region synchronization evaluation results.

[0151] First, the application performs a discrete Fourier transform (DFT) on the raw EEG data obtained in step 1, maps the time domain signal to the frequency domain, and obtains its spectrum representation. Specifically, assuming that the time series data of the jth channel of the i-th epoch is xi,j(n) (n=0,1,...,N-1), its spectrum Xi,j(k) can be calculated by DFT:

[0152] Xi,j(k) = Σ(n=0 to N-1) xi,j(n) * exp(-j2πkn / N), k=0,1,...,N-1.

[0153] Where N is the length of the time series, which is usually taken as an integer power of 2 to accelerate FFT calculation. For the complex frequency response Xi,j(k), this application focuses on its amplitude spectrum |Xi,j(k)| and phase spectrum angle(Xi,j(k)), which respectively characterize the intensity and phase characteristics of the channel at different frequency components. By splicing the spectrum of 64 channels, the overall frequency response diagram of the epoch can be obtained, that is:

[0154] Xi(k) = [|Xi,1(k)|*exp(jangle(Xi,1(k))), ..., |Xi,64(k)|*exp(jangle(Xi,64(k)))].

[0155] Among them, k=0 corresponds to the DC component, k=N / 2 corresponds to the Nyquist frequency (the highest resolvable frequency), and the rest are the fundamental wave and its multiple components.

[0156] Exemplarily, the obtaining of frequency domain feature information corresponding to the original EEG data includes: performing Fourier transform on the original EEG data to obtain a spectrum representation; generating a frequency response graph based on the spectrum representation; and extracting the frequency domain feature information based on the frequency response graph, wherein the frequency domain feature information includes at least a maximum response value and a frequency response flatness.

[0157] After obtaining the frequency response graph Xi(k), the present application can extract many useful frequency domain features, such as maximum response value, power spectrum distribution, band response, flatness, etc., to describe the frequency domain properties of the epoch:

[0158] 1. Maximum response value Amax,i = max(|Xi(k)|), k=0,1,...,N-1. Amax,i reflects the total energy or amplitude peak of the epoch and can be used to identify abnormal phenomena such as high-energy transients or mutations.

[0159] 2. Dominant frequency fmax,i = argmax(|Xi(k)|). fmax,i represents the dominant frequency component of the epoch and is an important reference for constructing a bandpass filter. If fmax,i falls within the known physiological frequency band (delta / theta / alpha / beta / gamma), it means that the epoch is mainly dominated by the corresponding EEG oscillation component.

[0160] 3. Frequency band response Br,i = Σ(fl <= |fk| <= fh) |Xi(k)|, r = delta / theta / ... According to the known physiological frequency band division, the present application can calculate the total response energy Br,i in each frequency band to evaluate the relative contribution of each oscillation mode.

[0161] 4. Power spectrum entropy Hs,i = -Σ(k=0 to N-1) p(k)logp(k), where p(k)=|Xi(k)|^2 / Σ|Xi(k)|^2 is the normalized power spectrum. Hs,i reflects the distribution of energy in different frequency components. The larger the value, the more uniform the energy distribution and the more complex the signal.

[0162] 5. Spectral flatness Fs,i = (exp(Hs,i) / N). The value of Fs,i is between 0 and 1. When Fs,i approaches 1, it indicates that the energy distribution in each frequency is uniform and the signal presents a flat frequency response characteristic; when Fs,i approaches 0, it indicates that the energy is highly concentrated in a few frequency components and the signal has a prominent frequency band characteristic. The present application can detect the imbalance of the distribution of frequency response in each frequency band based on Fs,i.

[0163] 6. Exponential flatness Es,i = Σ(k=0 to L-1) |Xi(k)|^a / Σ(k=L to N-1) |Xi(k)|^a, a>0. Es,i is used to characterize the imbalance of energy distribution between low-frequency and high-frequency components. When Es,i is far away from 1, it means that there is a significant difference in energy distribution between low-frequency and high-frequency components. By properly selecting L and a, this application can find abnormal energy distribution in a specific frequency band.

[0164] 7. Minimum description length MDL,i=-logP(Xi(k)|θ)+C(θ). The MDL criterion combines the likelihood term and the complexity penalty term, and can select a model θ with a simple structure and concentrated frequency response while maintaining a high likelihood. By minimizing MDL,i, the application can obtain the best model θ* for this epoch and its parameterized representation, which provides a starting point for frequency domain filtering and noise reduction.

[0165] The above frequency domain features comprehensively describe the frequency response characteristics of EEG signals, and provide a powerful tool for analyzing the frequency band distribution of signals, evaluating the degree of frequency response imbalance, detecting abnormal frequency components, etc. Combined with time domain features (such as periodic parameters of steps, distortion measures of steps, etc.), this application can more comprehensively grasp the behavior of EEG data in the time and frequency dimensions.

[0166] For example, for an epoch containing a typical alpha (8-13Hz) EEG component, if the present application finds that its band response Balpha is high and the dominant frequency fmax falls within the alpha band, while the power spectrum entropy Hs is low and the spectrum flatness Fs is low, it means that the energy of the epoch is highly concentrated in the narrow alpha band and can be regarded as an "ideal" alpha EEG signal; on the contrary, if Hs is high, Fs is close to 1, and Es is far from 1, it means that the energy distribution of the epoch is seriously unbalanced in each frequency band and may be contaminated by large noise, requiring targeted denoising. In addition, for epochs with abnormally large Amax or fmax deviating from the physiological frequency band, the present application will also mark them as data segments suspected of being affected by noise. By combining time-frequency domain features, the present application can efficiently identify suspicious distorted areas in EEG data.

[0167] At the same time, the differences in frequency response characteristics, data amplitude and periodic variation obtained from the steps extracted from different brain regions are compared to evaluate the degree of synchronization between different brain regions as auxiliary information for noise separation.

[0168] First, define the amplitude difference measure dA(i,j), period difference measure dT(i,j) and frequency response difference dF(i,j) between channels i and j: dA(i,j)=||Ai-Aj|| measures the amplitude distance between the two channel data in the time domain, where Ai and Aj are the data amplitude vectors of channels i and j respectively, and ||·|| represents the L2 norm; dT(i,j)=min(|Ti-Tj|,2π / (|Ti-Tj|)) is the normalized period difference, Ti and Tj are the dominant periods of the corresponding channels, and the minimum value of the absolute value difference |Ti-Tj| and the inverse of the period 2π / (|Ti-Tj|) is taken; dF(i,j)=||Fi-Fj|| is the L2 distance of the frequency responses of the two channels, and Fi and Fj are the normalized frequency response spectra obtained through the steps.

[0169] Next, the three difference metrics in amplitude, period and frequency domain are combined to define the synchronization index sync(i,j) of channels i and j:sync(i,j)=exp(-αdA(i,j)-βdT(i,j)-γ*dF(i,j)) where α, β and γ are adjustable weight parameters used to balance the relative contributions of the three difference metrics to synchronization. The value range of sync(i,j) is (0,1]. When dA(i,j), dT(i,j), and dF(i,j) are all close to 0, sync(i,j) will approach 1, indicating that the two channels are highly consistent in amplitude, period, and frequency domain, showing strong synchronization; on the contrary, when the three difference metrics are large, sync(i,j) will approach 0, indicating that the two channels have significant differences in time-frequency amplitude and have lost synchronization. This application calculates the synchronization index of all N(N-1) / 2 pairs of channels and sorts them from high to low, which can be regarded as a synchronization sequence: (i1,j1), (i2,j2), ..., (iM,jM) where M=N(N-1) / 2 is the total number of channel pairs, the top (ik,jk) represents a highly synchronized pair of brain regions, and the bottom (il,jl) may be affected by interference noise and lose synchronization.

[0170] Furthermore, the present application can threshold the synchronization sequence to obtain a synchronization matrix: syncMatrix(i,j)={1,if sync(i,j)≥τ 0,otherwise} where τ is the synchronization threshold, which can be pre-set or adaptively adjusted according to actual needs. This binary synchronization matrix intuitively reflects which brain region pairs are in a synchronized state (value 1) and which are in a desynchronized state (value 0) in the current time period. In the latter case, the desynchronized brain region pairs urgently need to eliminate the noise interference that causes desynchronization in order to restore the synchronization relationship they should have. Based on the synchronization matrix, the present application can construct an undirected synchronization graph model, treating brain regions as nodes and synchronization relationships as edge connections. In this way, with the help of graph theory algorithms and spectral clustering techniques, the community structure in the synchronization graph can be automatically discovered, that is, sub-clusters composed of several internally highly synchronized brain regions, which often correspond to neural network modules responsible for specific functions in the brain. Once a serious desynchronization phenomenon is detected in a sub-cluster, it means that the functional network is damaged and the normal neural oscillation is disturbed by noise. It is urgent to remove the noise component that causes desynchronization. This application can use the internal synchronization strength of each sub-cluster as a cluster-level feature and combine it with other indicators to guide the subsequent layered and block filtering processing strategy.

[0171] In addition to considering the global synchronization matrix and graph model, this application also focuses on the desynchronization of local brain regions and their time evolution. This application defines the instantaneous desynchronization degree of brain regions i and j at each time point t as: Dysync(i,j,t)=1-sync(i,j,t), if sync(i,j,t)<τ =0, otherwise, where sync(i,j,t) is the instantaneous synchronization of the i-th and j-th channels at time t, and Dysync(i,j,t) is assigned a non-zero value only when sync(i,j,t) is lower than the threshold τ. This indicator can accurately depict the occurrence, development and disappearance of the desynchronization phenomenon. For example, assuming that at t=10, Dysync(1,5,10) suddenly rises to 0.8, it means that channels 1 and 5 change from a synchronized state to a desynchronized state, which may be caused by a sudden increase in noise; if Dysync(1,5,t) continues to maintain a high level in the following time period, such as 0.8~0.9, it means that the desynchronized state is still continuing and a long period of filtering is required to restore synchronization; on the contrary, if Dysync(1,5,t) decays to 0 quickly, it can be judged as a short-term transient noise interference, and only local filtering of the data in this time period is required.

[0172] By accumulating Dysync(i,j,t) over time, the present application can also obtain the desynchronization duration Tdysync(i,j) and the desynchronization intensity integral Idysync(i,j), which respectively reflect the duration and cumulative intensity of the desynchronization event: Tdysync(i,j)=∫Dysync(i,j,t)dt Idysync(i,j)=∫Dysync(i,j,t)^2dt

[0173] These two indicators can help the present application identify pairs of brain regions that are subject to long-term or severe noise pollution. If the Tdysync of a pair of brain regions is large, it means that they are in a state of desynchronization for a long time, which may be due to the lasting influence of noise; if Idysync is large, it means that the desynchronization intensity has accumulated a lot of energy, and the noise interference has seriously distorted the normal coupling between brain regions. More filtering resources need to be invested to restore their synchronization state. For example, if the Tdysync (3, 17) and Idysync (3, 17) values ​​of channels 3 and 17 are both large, it means that the two brain regions have maintained a serious state of desynchronization for a long time, which may be due to some persistent noise source that has caused continuous interference to them, such as environmental electromagnetic radiation, physiological clutter, etc. In this case, the present application needs to perform targeted deep filtering on channels 3 and 17 in the step, and invest more computing power until the noise source is successfully blocked and the normal synchronous coupling of the two brain regions is restored.

[0174] Conversely, if the present application observes that both Tdysync(5,8) and Idysync(5,8) are maintained at a low level, and Dysync(5,8,t) has no obvious peak in the entire time course, it can be judged that channels 5 and 8 always maintain good synchronization, and there is no need to perform special filtering on them. This classification and targeted processing of out-of-synchronization brain area pairs can maximize the retention of effective EEG signals while efficiently eliminating noise interference, thereby improving the overall denoising effect.

[0175] Step S108, based on a preset causal reasoning framework, a target filtering algorithm is constructed according to the noise information and the degree of distortion, and interference noise corresponding to the original EEG data at each level is filtered out according to the target filtering algorithm.

[0176] Specifically, the method includes: (1) constructing a hierarchical causal graph model. First, the hierarchical division results {L1, L2, ..., LL} obtained in the steps are transformed into a directed hierarchical graph G = {V, E}. Where V is the set of nodes, and each node vi ∈ V corresponds to a brain region (or feature dimension) belonging to a specific level Lk. The edge set E contains two types of directed edges: (a) Intra-level edges: If nodes vi and vj belong to the same level Lk and there is a direct causal relationship vi → vj in the causal model of the steps, then an intra-level edge (vi, vj) is added to G; (b) Inter-level edges: If vi ∈ Lk and vj ∈ Ll (k < l) and there is a causal path starting from vi and terminating at vj in the causal model of the steps, then an inter-level edge (vi, vj) is added to G. The hierarchical causal graph G constructed in this way contains the causal relationships within each level and the causal transmission links between levels.

[0177] (2) Implanting a noise model on the hierarchical causal graph G. In this application, the hierarchical noise impact models NMi = {(Lk, NSk)|NSk = {noise type: (n1, p1), (n2, p2), ... (nP, pP)}} obtained in the steps are embedded into G, and each noise type np (1 ≤ p ≤ P) is simulated as a special noise node source_p. Specifically, for each noise type (np, pp) in the noise source distribution NSk of each Lk, this application adds a noise node source_p to G and connects a directed edge (source_p, vi) to all nodes vi in Lk, and sets the edge weight value of (source_p, vi) to pp. In this way, this application successfully embeds the noise model into the hierarchical causal graph G. The noise node source_p can be regarded as a virtual source injecting noise into the Lk layer, and the edge weight pp reflects the influence degree of the noise in this layer.

[0178] (3) Estimate variational parameters based on noise estimation and distortion measure. This application combines the results of steps and estimates the noise level η(vi) and distortion degree ρ(vi) of each node vi on the causal graph G. Where η(vi) is equal to the weighted sum of various noise influences on the level Lk where vi is located in the step noise influence model NMi: η(vi) = Σ(np, pp)∈NSk pp; and ρ(vi) is directly taken from the distortion measure Distortion(Lk) of the level Lk where vi is located in the step. Next, based on η(vi) and ρ(vi), this application constructs a variational optimization objective function: J(θ) = λΣ(vi∈V) ρ(vi)||x(vi)- f(θ, x(pa(vi)))||^2 + (1-λ)Σ(vi∈V) η(vi)||y(vi) - g(θ, y(pa(vi)))||^2 where x(vi) and y(vi) represent the original signal and noise component corresponding to node vi, respectively; f(θ, x(pa(vi))) and g(θ, y(pa(vi))) are inputs based on the parent node pa(vi), which predict the pure signal and noise component at vi, respectively, determined by the parameter θ; ||·|| represents the L2 norm; λ∈[0,1] is the adjustment parameter. The optimization objective function aims to minimize two parts: the first part is the reconstruction error weighted by the degree of distortion of node vi, which reflects the difference between the predicted value and the true pure signal x(vi) under the network weight θ; the second part is the noise residual weighted by the noise level of node vi, which reflects the difference between the predicted value and the true noise component y(vi). By optimizing the objective function J(θ), the present application can automatically learn the appropriate parameter θ*, which can effectively capture and remove the noise component y(vi) while reconstructing the original signal x(vi) as accurately as possible.

[0179] (4) Establish a local graph model based on synchronization analysis. Based on the synchronization evaluation results in the steps, the present application can automatically identify local graph structures composed of highly synchronized brain regions, which correspond to neural network modules responsible for different functions in the brain. Specifically, the present application first constructs a weight graph T=(V, E') in G, where the node set V is the same as G, both corresponding to brain regions or feature dimensions; the difference is that the weight of each undirected edge (vi, vj) in the edge set E' is equal to 1-sync(i,j), reflecting the degree of desynchronization of vi and vj. In this way, the originally highly synchronized brain region pair (vi, vj) will have a lower edge weight in T; while the desynchronized brain region pair will have a higher weight.

[0180] Next, the present application uses the spectral clustering algorithm to perform community discovery on T and divide different local network structures. In short, spectral clustering divides node clusters with high cohesion and low outgroup properties based on the eigenvector matrix of T. The nodes contained in these node clusters correspond to highly synchronized brain regions, and they have close network connections, which can quickly form stable signal transmission pathways. However, within the same cluster, if there are obvious desynchronization edges (i.e., high-weight edges), it means that the pathway is contaminated by some suspicious noise and requires careful denoising. To this end, the present application records each local network subgraph obtained by spectral clustering as {S1, S2, ..., SM}, and constructs M local optimization subproblems based on S1~SM, and calculates the optimization objective function Jm(θ) inside each subgraph Sm, and then solves the optimal weight parameter θm*. In this way, even if there is a serious desynchronization phenomenon in the local network, it can be accurately discovered and corrected, thereby retaining the neural oscillation pattern and signal quality inside the synchronous network to the greatest extent.

[0181] (5) Layered and block filtering decoupling. After the above four steps, this application is ready and can start to perform targeted layered and block filtering on different levels and different network modules to accurately remove various noise components. The specific filtering algorithm is as follows:

[0182] First, the present application divides the entire hierarchical causal graph G into L subgraphs {G_1, G_2, ..., G_L}, where G_k corresponds to the kth level L_k, including all nodes in L_k and their intra-layer and inter-layer edge connections. For each G_k, the present application combines it with the previously obtained local synchronization network {S_1, S_2, ..., S_M} to obtain a series of combined subgraphs {C_1, C_2, ..., C_N}. Each C_i contains some nodes in G_k and the overlapping local synchronization network S_j, which can simultaneously reflect the hierarchical causal relationship and synchronization properties of the nodes inside the subgraph.

[0183] Next, this application constructs the following hierarchical filtering model for each spliced ​​subgraph C_i:

[0184] minimizeΘ L(Θ) = Σ(v∈C_i) [λ_v||x(v) - f(Θ, x(pa(v)))||^2 + (1-λ_v)||y(v) - g(Θ, y(pa(v)))||^2] st λ_v = (1-η(v))ρ(v) / [η(v)(1-ρ(v)) + (1-η(v))ρ(v)].

[0185] Where Θ is the filter network parameter to be optimized; x(v) and y(v) are the original signal and noise component at node v, respectively; f(Θ, x(pa(v))) and g(Θ, y(pa(v))) are the predictions of the clean signal and noise component at v, respectively, based on the input of the parent node pa(v); λ_v is the weight coefficient dynamically calculated for each node v, which is determined by the noise level η(v) and distortion level ρ(v) of v. This multi-objective optimization problem aims to minimize the reconstruction error and noise residual of the clean signal at the same time, where the two loss terms of different nodes are dynamically weighted in the ratio of λ_v:(1-λ_v). Specifically, for a node v with a high noise level η(v) and a low distortion level ρ(v), its λ_v will approach 1, and the optimization goal will focus on removing the noise component y(v) of the node; conversely, for a node with a low η(v) and a high ρ(v), the optimization goal will focus on retaining the original signal x(v). Through this adaptive weighting method, the present application can apply filtering of different strengths to different nodes within the same spliced ​​subgraph C_i to achieve refined block denoising.

[0186] When optimizing the hierarchical filtering model, this application uses the alternating direction multiplier algorithm (ADMM) to decompose the original problem into multiple sub-problems that are easier to solve, and use distributed computing resources to achieve efficient parallel solutions. In short, ADMM introduces backup variables and Lagrange multipliers to split the original joint optimization problem into the following form:

[0187] x^(k+1) = argmin_x L_ρ(x, y^k, λ^k) y^(k+1) = argmin_y L_η(x^(k+1),y, λ^k) λ^(k+1) = λ^k + γ(x^(k+1) - y^(k+1)).

[0188] Among them, L_ρ and L_η correspond to the signal reconstruction loss and noise residual loss respectively, x and y are backup variables, λ is the Lagrange multiplier, and γ>0 is the algorithm step parameter. Based on the strategy of data parallelism and model parallelism, this application solves the above ADMM sub-problems in parallel on a distributed computing cluster, and iteratively updates x, y and λ until the algorithm converges or meets the stopping condition. This distributed solution paradigm based on optimization theory can make full use of existing heterogeneous computing resources such as GPUs and TPUs, significantly improve filtering efficiency, and meet the stringent requirements of real-time denoising.

[0189] Through the above-mentioned layered block filtering, the present application finally obtains the decoded signal x^(v) and the decoded noise y^(v) of each spliced ​​sub-graph C_i. For the original signal x(v), the strategy of the present application is to retain the decoded signal x^(v); and for the decoded noise y^(v), the present application adopts the following different processing methods:

[0190] (a) If node v belongs to a local synchronous network S_j, this application will retain y^*(v) as potential asynchronous noise inside S_j for further processing;

[0191] (b) If node v does not belong to any local synchronization network, this application directly removes y^*(v) and regards it as external noise interference in the corresponding level L_k.

[0192] (6) Separation of asynchronous noise in local network. For the asynchronous noise y^*(v) inside the local synchronous network retained in step (5), this application will design an adaptive spatial filter based on subspace projection to orthogonalize it with the noise subspace, thereby separating the asynchronous components such as Gaussian white noise and physiological noise. The filter can work online without the need to train parameters in advance as in step (5), and has extremely high computational efficiency.

[0193] Specifically, for a given local synchronization network S_j, the present application can construct a synchronization correlation matrix R=Σ(i=1 to q) s_i s_i^T based on its topological structure, where q is the number of nodes in S_j, and s_i is the standardized synchronization pattern vector corresponding to the i-th node. This correlation matrix R is actually the covariance matrix of all synchronization pattern vectors, and its eigenvectors just correspond to the principal component decomposition inside S_j.

[0194] The present application performs eigendecomposition on R to obtain its eigenvalues ​​lambda_1 ≥ lambda_2 ≥ ... ≥ lambda_q and the corresponding eigenvectors {u_1, u_2, ..., u_q}. Since the desired synchronization mode should occupy the dominant characteristic component of R, the present application first determines a threshold r (0≤r≤q), and regards the first r eigenvectors {u_1, u_2, ..., u_r} as the basis of the synchronization subspace, and the other eigenvectors {u_(r+1), ..., u_q} as the basis of the non-synchronization (noise) subspace.

[0195] Next, for each data x(v) belonging to S_j of node v, this application first projects it into the synchronization subspace:

[0196] x_s(v) = Σ(i=1 to r) (u_i^T * x(v))u_i.

[0197] This step can retain the components in x(v) that are highly correlated with the expected synchronization pattern. Then, this application uses the residual component of x(v) in the synchronization subspace as the asynchronous noise y(v):

[0198] y(v) = x(v) - x_s(v).

[0199] Since x_s(v) and y(v) are orthogonal to each other, the present application can easily separate the desired synchronous component x_s(v) and the asynchronous noise y(v) from the original observation data x(v). This adaptive subspace projection method can make full use of the prior knowledge of the synchronous network, and automatically identify and remove any form of asynchronous noise disturbance, including Gaussian white noise, electromyographic noise, electrooculographic noise, etc., while retaining the key synchronous components.

[0200] After step (5) of layered and block filtering and step (6) of local network asynchronous noise separation, the present application has effectively removed most of the layer noise and asynchronous noise. However, the residual wavelet noise and short transient noise may still have a certain impact on the signal quality, so the present application needs to further process them in detail, which is a "no residual processing" finishing link in the entire noise elimination process.

[0201] Wavelet noise refers to the high-frequency and low-frequency disturbance components existing in all scales and directions of wavelet decomposition. In order to accurately remove these residual noises while retaining the key wavelet signal energy, this application specifically designs an adaptive wavelet filter based on improved empirical mode decomposition (ICEMD).

[0202] The working principle of the filter is as follows: first, the data processed in steps (5) and (6) are subjected to wavelet decomposition to obtain multiple wavelet components W1, W2, ..., WM of different scales and directions. For each wavelet component Wi, the filter inputs it into the ICEMD module for decomposition, adaptively constructs the local characteristic curve of Wi, and extracts the intrinsic component function {IMF1, IMF2, ..., IMFn} and the residual r of Wi from it. Among them, each IMFj corresponds to a different oscillation mode in Wi, and the residual r often corresponds to the noise component.

[0203] Next, the filter will automatically analyze the statistical characteristics of each IMFj, such as amplitude and frequency, and compare them with the parameter range in the pre-modeled noise pattern library to determine whether the IMFj is noise. For IMFj that is confirmed to be noise, the filter will directly remove it; for non-noise IMFj, it will be retained. After a round of iterative decomposition and noise identification, this application will obtain the filtered version Wi' of the wavelet component Wi. Wi' = Σ(j belongs to non-noise IMF) IMFj.

[0204] Repeat the above steps, perform independent ICEMD decomposition and noise removal on each wavelet component Wi, and finally the present application will obtain the filtering results {W1', W2', ..., WM'} of all wavelet components. Perform wavelet reconstruction on these filtered wavelet components, and the present application can obtain the EEG signal x^w(t) with wavelet noise removed.

[0205] In addition to wavelet noise, this application also needs to pay special attention to short-term transient noise, such as sudden Sprick noise, short-term myoelectric noise, etc. This type of noise has a short duration and concentrated energy, but its impact may be very severe and will seriously damage the local signal quality. In order to efficiently identify and separate these transient noises, this application uses online dictionary learning algorithms such as K-SVD to automatically construct an over-complete dictionary set D that is highly correlated with EEG signals and noise patterns from x^w(t).

[0206] The dictionary D contains a large number of atomic molecules d1, d2, ..., dp, each di is a normalized codeword that can efficiently represent the local pattern in x^w(t). Through sparse decomposition, the present application can represent x^w(t) on the dictionary D as a linear superposition of a small number of non-zero coefficients: x^w(t) ≈ Σ αi * di, ||α||0 << p. Where α is a sparse coefficient vector, and ||α||0 represents the number of non-zero elements in α. The key is that since the dictionary D contains the atomic molecules of the EEG signal and the noise pattern at the same time, the signal component and the noise component must be projected on different dictionary atoms, so that the signal and noise parts can be automatically separated by analyzing the index position of the non-zero coefficient α.

[0207] Specifically, this application adopts a smoothing optimization framework to model the sparse decomposition of x^w(t) on the dictionary D as the following optimization problem: min ||α||0 + λ||x^w(t) - Dα||2^2 + γ∑ψ(αi). The first term ||α||0 is the l0 norm of the coefficient vector α, reflecting its sparsity; the second term is the signal fidelity term, which ensures that the decomposition result is as close to the original signal as possible; the third term is the smoothing term, which promotes the non-zero components of the coefficient vector α to maintain local continuity by designing a suitable smooth penalty function ψ (such as the total variation function). By optimizing and solving this problem, this application will obtain the optimal sparse decomposition α* and its corresponding signal component x^s and noise component x^n, that is:

[0208] x^w(t) = x^s + x^n = Σ(i belongs to the signal index) αi di + Σ(j belongs to the noise index) αjdj. The signal component x^s retains all the useful signal energy in x^w(t), while the noise component x^n includes various forms of transient noise components, such as sudden Sprick noise.

[0209] Finally, the present application combines and reconstructs the signal components x^s and x^w(t) after wavelet filtering and sparse decomposition to obtain the final pure EEG data x^*(t). x^*(t) = x^s + x^w(t). Through the above-mentioned "no residual processing" of time-frequency domain filtering and reconstruction, the present application effectively avoids further contamination of the pure signal x^(t) by various types of residual noise, while retaining the detailed features and local patterns in x^(t).

[0210] In order to implement the EEG data denoising method based on the hierarchical causal graph model corresponding to the above method embodiment, so as to achieve the corresponding functions and technical effects. Figure 2 , Figure 2 The structural block diagram of an EEG data denoising device 200 provided in an embodiment of the present application is shown. For the convenience of description, only the parts related to the present embodiment are shown. The EEG data denoising device 200 provided in an embodiment of the present application includes:

[0211] A data acquisition module 201 is used to acquire the original EEG data to be denoised and the time series corresponding to the original EEG data;

[0212] A data segmentation module 202 is used to segment the original EEG data according to the time series to obtain a plurality of segmented data and a starting point feature vector corresponding to each segmented data;

[0213] A model building module 203 is used to build a causal graph model based on a causal discovery algorithm and according to the starting point feature vector corresponding to each segmented data;

[0214] A change acquisition module 204, used for acquiring periodic change information corresponding to the original EEG data according to the causal graph model;

[0215] A hierarchical division module 205 is used to divide the causal graph model into multiple hierarchies; each hierarchy is used to reflect a different causal mechanism;

[0216] A noise acquisition module 206, configured to acquire noise information corresponding to each of the levels according to the causal graph model and preset prior knowledge;

[0217] A consistency acquisition module 207, configured to acquire data consistency information of the plurality of levels, and acquire the degree of distortion of each level according to the periodic variation information and the data consistency information;

[0218] The denoising completion module 208 is used to construct a target filtering algorithm based on a preset causal reasoning framework, according to the noise information and the degree of distortion, and filter out the interference noise corresponding to the original EEG data at each level according to the target filtering algorithm.

[0219] The above-mentioned EEG data denoising device 200 can implement the EEG data denoising method based on the hierarchical causal graph model of the above-mentioned method embodiment. The options in the above-mentioned method embodiment are also applicable to this embodiment and will not be described in detail here. The rest of the contents of the embodiment of the present application can refer to the contents of the above-mentioned method embodiment, and will not be repeated in this embodiment.

[0220] Figure 3 This is a schematic diagram of the structure of a computer device provided in one embodiment of the present application. Figure 3 As shown, the computer device 3 of this embodiment includes: at least one processor 30 ( Figure 3 Only one is shown in the figure), a memory 31 and a computer program 32 stored in the memory 31 and executable on the at least one processor 30, wherein the processor 30 implements the steps of any of the above method embodiments when executing the computer program 32.

[0221] The computer device 3 may be a computing device such as a smart phone, a tablet computer, a desktop computer, a cloud server, etc. The computer device may include but is not limited to a processor 30 and a memory 31. Those skilled in the art will understand that Figure 3 It is only an example of computer device 3 and does not constitute a limitation on computer device 3. It may include more or fewer components than shown in the figure, or a combination of certain components, or different components, for example, it may also include input and output devices, network access devices, etc.

[0222] The processor 30 may be a central processing unit (CPU), or other general-purpose processors, digital signal processors (DSP), application-specific integrated circuits (ASIC), field-programmable gate arrays (FPGA) or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. A general-purpose processor may be a microprocessor or any conventional processor, etc.

[0223] In some embodiments, the memory 31 may be an internal storage unit of the computer device 3, such as a hard disk or memory of the computer device 3. In other embodiments, the memory 31 may also be an external storage device of the computer device 3, such as a plug-in hard disk, a smart media card (SMC), a secure digital (SD) card, a flash card, etc. equipped on the computer device 3. Further, the memory 31 may also include both an internal storage unit and an external storage device of the computer device 3. The memory 31 is used to store an operating system, an application program, a boot loader (BootLoader), data, and other programs, such as the program code of the computer program. The memory 31 may also be used to temporarily store data that has been output or is to be output.

[0224] In addition, an embodiment of the present application further provides a computer-readable storage medium, wherein the computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, the steps in any of the above method embodiments are implemented.

[0225] An embodiment of the present application provides a computer program product. When the computer program product is run on a computer device, the computer device implements the steps in the above-mentioned method embodiments when executing the computer device.

[0226] In several embodiments provided in the present application, it is understood that each box in the flow chart or block diagram can represent a module, a program segment or a part of a code, and the module, a program segment or a part of a code contains one or more executable instructions for realizing the specified logical function. It should also be noted that in some alternative implementations, the functions marked in the box can also occur in a different order from that marked in the accompanying drawings. For example, two consecutive boxes can actually be executed substantially in parallel, and they can sometimes be executed in the opposite order, depending on the functions involved.

[0227] If the functions are implemented in the form of software function modules and sold or used as independent products, they can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present application can be essentially or partly embodied in the form of a software product that contributes to the prior art. The computer software product is stored in a storage medium and includes several instructions for a computer device to execute all or part of the steps of the methods described in the various embodiments of the present application. The aforementioned storage media include: U disk, mobile hard disk, read-only memory (ROM, Read-Only Memory), random access memory (RAM, Random Access Memory), disk or optical disk, and other media that can store program codes.

[0228] The specific embodiments described above further describe the purpose, technical solutions and beneficial effects of the present application in detail. It should be understood that the above description is only a specific embodiment of the present application and is not intended to limit the scope of protection of the present application. It is particularly pointed out that for those skilled in the art, any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present application should be included in the scope of protection of the present application.

Claims

1. A method for denoising EEG data based on a hierarchical causal graph model, characterized in that: include: Acquire the original EEG data to be denoised and the time series corresponding to the original EEG data; Segmenting the original EEG data according to the time series to obtain a plurality of segmented data and a starting point feature vector corresponding to each segmented data; Based on the causal discovery algorithm, a causal graph model is constructed according to the starting point feature vector corresponding to each of the segmented data; Acquiring periodic variation information corresponding to the original EEG data according to the causal graph model; Dividing the causal graph model to obtain multiple levels; each level is used to reflect a different causal mechanism; According to the causal graph model and preset prior knowledge, obtaining noise information corresponding to each of the levels; Acquire data consistency information of the plurality of levels, and acquire the degree of distortion of each level according to the periodic variation information and the data consistency information; Based on a preset causal reasoning framework, a target filtering algorithm is constructed according to the noise information and the degree of distortion, and interference noise corresponding to the original EEG data at each level is filtered out according to the target filtering algorithm.

2. The method according to claim 1, characterized in that The causal discovery algorithm is based on which a causal graph model is constructed according to the starting point feature vector corresponding to each segmented data, including: Based on the skeleton recovery algorithm, an undirected skeleton graph is obtained in each of the starting point feature vectors; Inputting the undirected skeleton graph into a graph neural network encoder, and outputting a representation vector corresponding to the undirected skeleton graph; Sampling each of the representation vectors according to a Monte Carlo graph sampling strategy to obtain a plurality of candidate directed acyclic graph samples; Obtaining a candidate graph score corresponding to each candidate directed acyclic graph sample, and determining a target sample from the plurality of candidate directed acyclic graph samples according to the candidate graph score; The target sample is input into a preset graph neural network decoder, and the causal graph model is output.

3. The method according to claim 2, characterized in that The step of obtaining an undirected skeleton graph from each starting point feature vector includes: Constructing a completely undirected graph according to a plurality of starting point feature vectors; A conditional independence test is performed on the completely undirected graph, and edges that do not meet a preset separation criterion are deleted from the completely undirected graph to obtain the undirected skeleton graph.

4. The method according to claim 1, characterized in that: The step of obtaining periodic variation information corresponding to the original EEG data according to the causal graph model includes: Identifying a periodic embedding vector corresponding to each of the causal graph models according to a periodicity detection algorithm based on a graph convolutional neural network; Decoding is performed according to the periodic embedding vector to obtain a periodic parameter; The periodic variation information is generated according to the periodic parameter.

5. The method according to claim 1, characterized in that The causal graph model is divided to obtain multiple levels, including: Topologically sorting the causal graph model to obtain an ordered node sequence; Obtaining node path information corresponding to the ordered node sequence; The ordered node sequence is divided into a plurality of levels according to the node path information.

6. The method according to claim 1, characterized in that The obtaining of noise information corresponding to each of the levels includes: Obtaining connection topology attributes corresponding to each of the layers; Acquire the noise propagation path and impact range according to the connection topology attribute; The noise information is generated according to the noise propagation path and the influence range.

7. The method according to claim 1, characterized in that Before constructing the target filtering algorithm according to the noise information and the degree of distortion, the method further includes: Acquire frequency domain feature information corresponding to the original EEG data; According to the frequency domain feature information and the periodic variation information, the brain region synchronization evaluation result is obtained, so as to construct a target filtering algorithm according to the noise information, the degree of distortion, the target filtering algorithm and the brain region synchronization evaluation result.

8. The method according to claim 7, characterized in that The obtaining of frequency domain feature information corresponding to the original EEG data includes: Perform Fourier transform on the original EEG data to obtain the spectrum representation; generating a frequency response graph based on the frequency spectrum representation; The frequency domain characteristic information is extracted according to the frequency response diagram, and the frequency domain characteristic information at least includes a maximum response value and a frequency response flatness.

9. An EEG data denoising device based on a hierarchical causal graph model, characterized in that: include: A data acquisition module, used to acquire the original EEG data to be denoised and the time series corresponding to the original EEG data; A data segmentation module, used to segment the original EEG data according to the time series, and obtain a plurality of segmented data and a starting point feature vector corresponding to each segmented data; A model building module, used to build a causal graph model based on a causal discovery algorithm and according to the starting point feature vector corresponding to each segmented data; A change acquisition module, used for acquiring periodic change information corresponding to the original EEG data according to the causal graph model; A hierarchical division module, used to divide the causal graph model into multiple hierarchies; each hierarchy is used to reflect a different causal mechanism; A noise acquisition module, used for acquiring noise information corresponding to each of the levels according to the causal graph model and preset prior knowledge; A consistency acquisition module, used to acquire data consistency information of the multiple levels, and acquire the degree of distortion of each level according to the periodic variation information and the data consistency information; The denoising completion module is used to construct a target filtering algorithm based on a preset causal reasoning framework according to the noise information and the degree of distortion, and filter out the interference noise corresponding to the original EEG data at each level according to the target filtering algorithm.

10. A computer device, characterized in that: The method comprises a memory and a processor; the memory is used to store a computer program; the processor is used to execute the computer program and implement the method according to any one of claims 1 to 8 when executing the computer program.

Images