Second-order blind identification multi-artifact suppression optimization method and system based on optical pump magnetometer

Abstract

The invention discloses a second-order blind identification multi-artifact suppression optimization method and system based on an optical pump magnetometer, and belongs to the technical field of biomedical crossing field signal processing. The method comprises the following steps: acquiring multi-channel OPM-MEG data, wherein the multi-channel OPM-MEG data comprises brain signal data and artifact data; on the basis of the collected data, combined approximate diagonalization of a single-delay covariance matrix is carried out, a sensitivity index is constructed, the contribution degree of a specific delay point to separated data is evaluated, and a data-driven parameter adjustment strategy is generated; based on artifact priori knowledge, respectively constructing artifact templates for low-specificity artifacts and high-specificity artifacts to form a target function; and carrying out SOBI source separation on the acquired multi-channel OPM-MEG data, and optimizing the separation process based on the generated parameter adjustment strategy and the target function to realize artifact suppression. According to the method, the problem of parameter adaptive selection of the SOBI in different artifact scenes can be solved.

Description

Technical Field

[0001] The present invention belongs to the technical field of signal processing in the biomedical cross - field, and particularly relates to an optimized method and system for second - order blind identification and multi - artifact suppression based on an optically pumped magnetometer, which is applicable to the simultaneous suppression of multiple artifacts in magnetoencephalography (MEG). Background Art

[0002] Magnetoencephalography (MEG) is a non - invasive brain functional imaging method with high spatio - temporal resolution, and has made remarkable progress in the fields of brain function, cognitive research, and clinical applications of epilepsy. In recent years, the application of optically pumped magnetometers (OPMs) has overcome the dependence of traditional superconducting quantum interference devices (SQUIDs) on cryogenic environments, can be placed closer to the scalp, improve the signal - to - noise ratio of measured signals, and potentially enhance the spatial resolution of source reconstruction. However, MEG signals are weak and vulnerable to various artifacts, such as low - frequency interference caused by physiological artifacts, metal artifacts, and motion artifacts. Removing these interferences in a multi - artifact scenario is a major challenge in MEG data processing. Currently, there are various methods for removing artifacts from MEG data. Among them, the blind source separation (BSS) algorithm can simultaneously separate multiple artifact components from complex mixed signals without prior knowledge of the mixing process or source signal characteristics, significantly reducing the need to separately set reference channels for each artifact. In addition, it also shows strong adaptability when facing complex artifacts, especially suitable for complex scenarios with co - existing multiple artifacts, so it has become the mainstream choice for MEG artifact removal.

[0003] BSS algorithms include FastICA, Infomax, AMUSE, JADE, and SOBI, etc. Among them, SOBI has received attention due to its excellent performance in the separation and identification of complex neural activity patterns.

[0004] However, the performance of SOBI highly depends on the selection of the time delay set τ, which directly determines the construction of the optimization objective in joint approximate diagonalization and thus affects the separation effect. Currently, common methods include repeatedly trying different time delay combinations or selecting a relatively long and complete time delay set to ensure that the frequency components of the source signals can fall within the frequency bands related to the time delay values. The former lacks an objective evaluation criterion, and the selection process relies on the subjective judgment of the researcher, which is time-consuming and laborious. Although the latter is conservative, the large number of introduced time delay parameters will increase the convergence difficulty of joint diagonalization and reduce the independence of the decomposed components. In addition, adjacent time delay points provide highly redundant information. For each additional time delay point, new noise components may be introduced, and their cumulative effect will obscure the characteristics of the true signal components. In practical OPM-MEG applications, the variation of artifact data characteristics is common and inevitable. Different experimental conditions (such as noise level, environmental temperature, device performance) and the physiological states of participants (such as fatigue, mood swings, task paradigms) will all affect the signals. The basic physiological artifacts of participants at different times will also show differences due to various factors, and the differences are more significant among different individuals. In addition, there are complex artifact sources such as metal implants, which will be modulated by types and other factors, making the situation more complex and variable. Therefore, the fixed and conservative parameter selection method lacks adaptability and cannot meet the requirements of variable scenarios. Summary of the Invention

[0005] To solve the above technical problems, the present invention provides an optimization method and system for second-order blind identification and multi-artifact suppression based on an optically pumped magnetometer, which can dynamically adjust the time delay parameters to adapt to different scenarios and artifact characteristics. By fully utilizing the data characteristics of the source decomposed by a single time delay point, time delay redundancy removal and partitioning are performed, and a guiding template is constructed in combination with the prior knowledge of artifacts. A dynamic feedback mechanism is established, and parameter adjustment is carried out through dynamic evaluation of the separation effect. Under a specific setting of the number of delay points, the time delay combination that can best distinguish the target artifacts is found, and various known artifacts are removed to the greatest extent while retaining the neural signals. Simulation experiments and auditory evoked experiments in different artifact data scenarios in the OPM-MEG environment verify the effectiveness of the proposed method. Compared with the fixed time delay set, this method can flexibly adjust the time delay parameters when facing different degrees and types of artifact interferences, continuously improving the signal-to-noise ratio and source localization accuracy.

[0006] To achieve the above object, the present invention adopts the following technical solutions:

[0007] In a first aspect, the present invention provides an optimization method for second-order blind identification and multi-artifact suppression based on an optically pumped magnetometer, including the following steps:

[0008] Step 1, collect multi-channel OPM-MEG data, where the multi-channel OPM-MEG data includes brain signal data and artifact data;

[0009] Step 2: Based on the collected data, perform joint approximate diagonalization of the single-delay covariance matrix, construct sensitivity indicators, evaluate the contribution degree of specific delay points to the separated data, and generate a data-driven parameter adjustment strategy;

[0010] Step 3: Based on the prior knowledge of artifacts, construct artifact templates for low-specificity artifacts and high-specificity artifacts respectively to form an objective function;

[0011] Step 4: Perform SOBI source separation on the collected multi-channel OPM-MEG data, and optimize the separation process based on the generated parameter adjustment strategy and the objective function to achieve artifact suppression.

[0012] In a second aspect, the present invention provides an optimization system for second-order blind identification and multi-artifact suppression based on an optically pumped magnetometer, including: a data acquisition module, a strategy construction module, an objective function construction module, and an optimization module. Among them:

[0013] The data acquisition module is used to acquire multi-channel OPM-MEG data, and the multi-channel OPM-MEG data includes brain signal data and artifact data;

[0014] The strategy construction module is used to perform joint approximate diagonalization of the single-delay covariance matrix based on the collected data, construct sensitivity indicators, evaluate the contribution degree of specific delay points to the separated data, and generate a data-driven parameter adjustment strategy;

[0015] The objective function construction module is used to construct artifact templates for low-specificity artifacts and high-specificity artifacts respectively based on the prior knowledge of artifacts to form an objective function;

[0016] The optimization module is used to perform SOBI source separation on the collected multi-channel OPM-MEG data, and optimize the separation process based on the generated parameter adjustment strategy and the objective function to achieve artifact suppression.

[0017] In a third aspect, the present invention provides an electronic device, including: one or more processors; a memory for storing one or more programs; wherein, when the one or more programs are executed by the one or more processors, the one or more processors implement the foregoing method for second-order blind identification and multi-artifact suppression based on an optically pumped magnetometer.

[0018] In a fourth aspect, the present invention provides a computer-readable storage medium, on which executable instructions are stored, and when the instructions are executed by a processor, the processor can implement the foregoing method for second-order blind identification and multi-artifact suppression based on an optically pumped magnetometer.

[0019] The beneficial effects of the present invention are as follows:

[0020] The present invention proposes an optimized framework for multi-artifact suppression in OPM-MEG based on second-order blind identification, aiming to address the challenges in the adaptive selection of parameters in the SOBI algorithm. By integrating data-driven sensitivity analysis, prior-knowledge-guided template construction, and a dynamic feedback mechanism, a transformation from "blind" source separation to "informed" source separation is achieved;

[0021] The present invention defines a sensitivity index and a selection strategy for evaluating the quality of a single time-delay correlation matrix, and based on this, selects the delay points that reflect important scale information while eliminating the low-contribution delay points to accelerate the optimization process, improve the separation efficiency and the artifact suppression effect.

[0022] The present invention relies on the detailed modeling of various artifacts in magnetoencephalography experiments to form a set of standard template systems so that more types of artifacts can be incorporated in the future. Brief Description of the Drawings

[0023] Figure 1 is a framework diagram of the optimized method for multi-artifact suppression based on second-order blind identification using an optically pumped magnetometer in the present invention, where (a) is the data-driven and parameter adjustment part, (b) is the prior integration and artifact template part, and (c) is the parameter automatic optimization framework part;

[0024] Figure 2 is a general diagram of the analysis of multi-artifact high specificity and low specificity. Detailed Embodiment

[0025] The present invention will be further described below with reference to the drawings and embodiments.

[0026] The present invention improves the efficiency and accuracy of parameter selection and the artifact suppression effect by integrating data-driven sensitivity analysis, prior-knowledge-guided template construction, and a dynamic feedback mechanism. In particular, the method can adaptively adjust parameters according to different scenarios and corresponding prior knowledge, rather than relying on multiple attempts or using a large set of time delays, providing the possibility for rapid parameter adjustment according to scenarios in the future. The simulation experiment results show that with fewer optimized delay points, the method achieves better signal separation effect and source localization accuracy than classical parameters, while reducing the demand for computing resources.

[0027] Specifically, as Figure 1 shown, it includes the following steps:

[0028] Step 1, collect multi-channel OPM-MEG data, where the multi-channel OPM-MEG data includes brain signal data and artifact data;

[0029] Step 2: Based on the collected sensor data, perform joint approximate diagonalization of the single-delay covariance matrix, construct sensitivity indicators, evaluate the contribution degree of specific delay points to the separated data, and generate a data-driven parameter adjustment strategy, such as Figure 1 shown in (a) of

[0030] Step 3: Based on the prior knowledge of artifacts in the mixed data, construct artifact templates for low-specificity artifacts and high-specificity artifacts respectively, form an objective function for judging the accuracy of blind separation, such as Figure 1 shown in (b) of

[0031] Step 4: Perform SOBI source separation on the collected OPM-MEG data, optimize the separation process based on the generated parameter adjustment strategy and the constructed objective function, and achieve artifact suppression, such as Figure 1 shown in (c) of

[0032] In the above Step 1, multi-channel OPM-MEG data is collected. The multi-channel OPM-MEG data includes brain signal data and artifact data, including collecting brain signal data using a multi-channel OPM sensor array and collecting eye movement, heartbeat, and metal artifact data using a specific reference sensor. The data acquisition model is , where represents the data recorded by the OPM sensor, represents the number of OPM sensors, represents the number of sampling points, is the source activity signal, is the number of source activity signals, represents the linear mixing matrix, describing the mixing mapping from each source activity to the sensor space, represents random noise.

[0033] In the above Step 2, for delay points, the possible number of delay combinations is . In the case of high-sampling-rate data, the exponential growth of this combination number exceeds the processing capacity of the current computer and cannot be traversed for optimization. Therefore, based on the collected sensor data, perform joint approximate diagonalization of the single-delay covariance matrix, construct sensitivity indicators, evaluate the contribution degree of specific delay points to the separated data, and generate a data-driven parameter adjustment strategy to accelerate the calculation and improve the separation effect. It includes:

[0034] Step 2.1: Based on the collected sensor data , construct a single-delay covariance matrix . The joint approximate diagonalization operation is to find a common orthogonal transformation matrix to transform the covariance matrix into a diagonal matrix. The optimal transformation matrix where τ represents a single time delay, represents the sum of the squares of the non-diagonal elements of the matrix, is the transpose rank.

[0035] Step 2.2: By constructing a sensitivity index , quantify the diagonalization trend of the single time delay matrix and the imbalance of the source signal energy distribution , to evaluate the contribution of a specific single time delay to separating the current data, which is defined as follows:

[0036] (1)

[0037] (2)

[0038] (3)

[0039] where norm represents the normalization operation, max(X) represents the maximum element value in matrix X, is the matrix eigenvalue matrix, which is sorted in descending order by variance by default, represents the i-th eigenvalue of the eigenvalue matrix, q is the number of eigenvalues taken with a fixed quantity, the superscript T represents the transpose, m is the total number of eigenvalues, and its value is equal to the number of OPM sensors ;

[0040] Step 2.3: Calculate the sensitivity index for each possible single time delay , and perform smoothing processing using Gaussian filtering; adopt a multi-scale peak detection method to identify the peaks of the sensitivity index at different time scales, and form a peak group by taking the preset proportion of delay points around the peak; use the DBSCAN clustering algorithm to identify different peak groups and group them into intervals, and the set of intervals is denoted as , where K represents the total number of intervals;

[0041] Step 2.4: Based on the divided set of intervals, construct a data-driven parameter adjustment strategy. Considering the sensitivity difference of the source signal to different time delays, introduce a selection probability P for selecting different delay intervals; introduce an interval weight W to assign different weights to intervals containing different time scales; combine the interval selection probability P and the interval weight W, select multiple intervals and perform random sampling of delay points respectively, thereby forming an optimization constraint condition for parameter adaptive adjustment and forming a complete parameter adjustment strategy.

[0042] Construct the constraint conditions for strategy selection as follows:

[0043] (4)

[0044] (5)

[0045] (6)

[0046] (7)

[0047] in, represents the jth interval The number of discrete points contained in , q represents the power parameter of the nonlinear attenuation function, represents the set of available intervals when selecting the tth delay point, Indicates the interval selected when selecting the tth delay point, "\" indicates the difference operation of the set, represents the selected tth delay point, and U represents uniform distribution.

[0048] In the above step 3, in the actual MEG signal separation, completely separating all source signals is not only technically difficult (the components are complex and mostly non-stationary signals), but also not necessary (the specific meaning of each unknown component cannot be accurately identified); instead, the focus should be on the extraction of key signals or the suppression of specific artifacts. Therefore, based on the prior knowledge of artifacts in mixed data, artifact templates are constructed for low-specificity artifacts and high-specificity artifacts respectively to form an objective function for the accuracy judgment of blind separation, thereby helping to judge the accuracy of "blind" separation and improving the interpretability of decomposed components. Including:

[0049] Step 3.1: Artifacts are divided into low-specificity artifact signals with relatively fixed characterizable signal characteristics or propagation modes, and high-specificity artifact signals that need to rely on reference sensors due to individual differences or large influencing factors. The analysis of artifact signal characteristics is summarized in Figure 2 ;

[0050] Step 3.2: For low-specificity artifacts, corresponding time, frequency, and space-domain feature templates can be established for identification and suppression. For example, for blink artifacts, the time domain feature is a single obvious instantaneous jump and rapid fall, and the activated area is located in the frontal lobe. For heartbeat artifacts, the time domain feature is a combination of QRS complexes and T waves, and the activated area is located in the left temporal lobe or occipital lobe. Based on this, the corresponding prior template can be constructed, and the Pearson correlation coefficient can be used for identification and score calculation. Considering that the pre-constructed artifact template cannot fully adapt to each individual, adaptive migration is required, and judgment is made on data that is difficult to extract artifacts. See the schematic diagram. Figure 1In (b), the specific steps include: clustering MEG sensors using prior knowledge of artifacts, compressing the scale to improve the signal-to-noise ratio of specific artifacts, matching the decomposition results of local data with a prior template to identify the most similar components; if the similarity is above the threshold, update the template for subsequent global identification; if it is below the threshold, it indicates that the artifact feature is not significant and is not considered.

[0051] Step 3.3: For high-specificity artifacts, it is necessary to further analyze their features and rely on additional reference channel information. For example, metal artifacts are affected by the type, location, and surrounding environment of the implanted metal, and their features vary greatly among different subjects. The low-frequency power ratio of 0.5 - 8 Hz can be used in combination with reference channel information, and mutual information is used to extract deep correlations for identification and score calculation.

[0052] Step 3.4: Calculate the matching scores for the identified potential various artifact components:

[0053] (8)

[0054] (9)

[0055] (10)

[0056] (11)

[0057] Among them, and respectively represent the matching scores of low- and high-specificity artifacts, represents the Pearson correlation coefficient, x and y respectively represent the source signal of the blind separation component and the template signal, represents the fast Fourier transform, represents the absolute value operation, a and b are the corresponding spatial activation patterns, represents the mutual information index, represents a specific quantization index constructed based on artifact analysis, and respectively represent the i-th value of the source signal x of the blind separation component and the template signal y, and are respectively the means of the source signal x of the blind separation component and the template signal y, n is the number of sampling points, p(x,y) represents the joint probability distribution, and p(x) and p(y) are respectively the marginal probability distributions of the source signal x of the blind separation component and the template signal y, and are respectively specific values of the source signal x of the blind separation component and the template signal y;

[0058] Step 3.5: Maximizing the objective function for separating known artifacts, that is, the optimization problem is:

[0059] (12)

[0060] (13)

[0061] Among them, and respectively represent the matching scores of the e-th low-specificity artifact and the f-th high-specificity artifact under a set of k delay parameters. is the penalty coefficient; N is the number of finally removed components; N max is the upper threshold of the removed components. represents the total set of optional time delays. Introducing equation is designed to avoid signal loss caused by excessive removal, and at the same time prevent noise from not being effectively separated and affecting the independence of other components. Specifically, when the number of removed components exceeds the upper threshold N max , the score will be non-linearly suppressed to balance the actual requirements between artifact removal and neural signal integrity.

[0062] In step 4, the collected OPM-MEG data B is first whitened to obtain the whitened data matrix Z; calculate the sensitivity index of the data Z at different time delay points, and perform interval division according to the index values to determine the initial delay point set; based on this delay point set, use the SOBI method to perform blind separation on Z, and optimize and solve the orthogonal transformation matrix V through Jacobi iteration to obtain the separated source signal components and their corresponding spatial activation patterns; combine the objective function constructed in step 3 to perform matching recognition and quantitative evaluation of the artifact components on the separation result; combine the parameter adjustment strategy generated in step 2 to perform directed dynamic reallocation on the delay combination; establish a feedback mechanism through an intelligent optimization algorithm to adjust and optimize the selection of the delay combination to improve the separation effect and quality in a multi-artifact interference environment.

[0063] To evaluate the effectiveness of the method proposed in the present invention, a series of semi-simulation experimental tests based on measured data were constructed. Multiple metrics were used to evaluate the algorithm performance, including root mean square error (RMSE), signal-to-noise ratio (SNR), and localization error (LE). Denote as the simulated brain internal signal, as the external noise signal, as the simulated signal detected by the sensor, and the signal after noise reduction and reconstruction is defined as .

[0064] The root mean square error (RMSE) reflects the distortion degree of the signal, and its calculation formula is:

[0065] (14)

[0066] The signal-to-noise ratio (SNR) reflects the clarity of the signal, and its calculation formula is:

[0067] (15)

[0068] Considering that the ultimate goal of MEG signal analysis is source localization, the localization error (LE) is introduced as an evaluation index. LE is defined as the Euclidean distance between the true simulated dipole position and the dipole position obtained at the signal peak moment using the equivalent current dipole (ECD) method.

[0069] To simulate the artifacts in the OPM-MEG measurement environment, the following simulations were carried out in the present invention: The heartbeat artifact was simulated by placing a dipole source half a meter to the left of the heart, simulating the heartbeat frequency of a normal adult (50 - 80 beats per minute); the blink artifact was simulated by placing symmetric dipole sources in the eye area and randomly generating blink events (15 - 30 times per minute); referring to the existing literature research, the metal artifact was simulated by a method based on measured data to simulate the influence of a metal tooth correction device. By measuring the resting state data of multiple batches of subjects with / without metal orthodontic brackets and , the spatial activation pattern of the main interference channels was analyzed to obtain the main interference component O, and the artifact propagation matrix was calculated, and the selected artifact components were projected onto all sensor channels of the artifact-free resting state using linear regression, that is to simulate the metal artifact interference process. A random dipole in the superior temporal gyrus region was selected as the simulated internal brain source for the neural source design. During the simulation experiment, different resting state data (including with and without metal artifacts), heartbeat artifacts, blink artifacts, and random noise were randomly combined, and a total of 15 different simulation scenario data (average SNR was -19.8 dB) were created for the following tests:

[0070] Five simulation scenario data were taken to test the effectiveness of the time delay criterion. The five groups of simulation data were optimized under 2, 4, 6, 8, and 10 time delay parameters respectively. Each group of data was repeated five times, the number of iterations was limited to three times, and the same component automatic rejection criterion was used. By using and not using the time delay partition criterion, the differences in the results were compared. The results are shown in Table 1, indicating that using time delay partition screening can find better values faster and more effectively in a short time compared to relying solely on the optimization algorithm, manifested as a higher signal-to-noise ratio (SNR) and a lower root mean square error (RMSE). These results show that time delay partition screening effectively improves the optimization speed and accuracy.

[0071] Table 1

[0072] Take the data of 1 simulation scenario and test the influence of the number of delay points in short-term iteration on the results. Under the condition of limiting three iterations, conduct multiple repeated experiments and adopt the same component rejection criteria. Use the adaptability score of the optimization objective as the comparison result, and the results are shown in Table 2.

[0073] Table 2

[0074] As can be seen from Table 2, different numbers of delay points result in different levels of optimization difficulty. As the number of delay points increases, the results tend to be stable. This is because the increase in information volume reduces the impact of internal variations at individual points on the overall result. However, this also makes the optimization process more complex and difficult to achieve the ideal result. On the contrary, a smaller number of delay points can obtain a smaller error value in a short time, but the variance increases, increasing the risk of difficult to find a better solution. Especially when the number of delay points is too small, the optimization process becomes extremely unstable and lacks robustness, making it difficult to meet the actual requirements.

[0075] Take the data of 9 simulation scenarios for comparative testing to evaluate the adaptive adjustment ability of the proposed method and the classical time-delay method in dealing with complex artifact interferences in different scenarios. Referring to the existing literature research, the currently used classical benchmark time-delay parameters are:

[0076] Time-delay set 1: {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 14, 16, 18, 20, 25, 35, 40, 45,50, 60, 70, 75, 80, 85, 90, 95, 100, 120, 140, 160, 180, 200, 240, 260, 280,300, 350, 400, 450, 500, 550, 600, 650, 700, 750, 800, 850, 900, 950, 1000,1100, 1200, 1300, 1400, 1500, 1600, 1700, 1800, 1900, 2000}

[0077] Time-delay set 2: {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 14, 16, 18, 20, 25, 30, 35, 40,45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100, 120, 140, 160, 180, 200,220, 240, 260, 280, 300}

[0078] Time delay set 3: {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50}

[0079] Table 3 shows the SNR and RMSE results at the sensor level after decomposition and component removal using different time delay parameters as inputs in each simulation scenario. The results show that the effect of the fixed time delay parameter method fluctuates greatly in different artifact scenarios, while the method of the present invention can achieve higher signal-to-noise ratio and lower error in most cases. The key to this stability and performance improvement lies in that the method of the present invention can dynamically adjust the time delay parameter according to data characteristics and prior guidance, so as to better adapt to different scenario data. However, since the optimization algorithm is carried out under short-term iteration and cannot fully explore all possible parameter spaces, there are performances in a few cases that fail to exceed the classical parameters. Increasing the number of iterations or introducing a more comprehensive parameter exploration mechanism may further improve the performance, but this will come at the cost of increased computing time.

[0080] Table 3

[0081] On the other hand, the present invention provides a second-order blind identification multi-artifact suppression optimization system based on an optically pumped magnetometer, and each module included therein can implement each step of the foregoing method. Specifically, it includes: a data acquisition module, a strategy construction module, an objective function construction module, and an optimization module. Among them:

[0082] The data acquisition module is used to acquire multi-channel OPM-MEG data, and the multi-channel OPM-MEG data includes brain signal data and artifact data;

[0083] The strategy construction module is used to perform joint approximate diagonalization of the single time delay covariance matrix based on the acquired data, construct a sensitivity index, evaluate the contribution degree of specific delay points to the separated data, and generate a data-driven parameter adjustment strategy, as Figure 1 shown in (a) below;

[0084] An objective function construction module, which is configured to construct artifact templates for low-specificity artifacts and high-specificity artifacts respectively based on artifact prior knowledge, and form an objective function, as shown in Figure 1 Figure (b) in

[0085] An optimization module, which is configured to perform SOBI source separation on the collected multi-channel OPM-MEG data, and optimize the separation process based on the generated parameter adjustment strategy and the objective function to achieve artifact suppression, as shown in Figure 1 Figure (c) in

[0086] In a third aspect, the present invention provides an electronic device, including: one or more processors; a memory for storing one or more programs; wherein, when the one or more programs are executed by the one or more processors, the one or more processors implement the foregoing second-order blind identification multi-artifact suppression optimization method based on an optically pumped magnetometer.

[0087] In a fourth aspect, the present invention provides a computer-readable storage medium, on which executable instructions are stored, and when the instructions are executed by a processor, the processor can implement the foregoing second-order blind identification multi-artifact suppression optimization method based on an optically pumped magnetometer.

[0088] The specific embodiments described above further elaborate on the purpose, technical solutions, and beneficial effects of the present invention. It should be understood that the above are only specific embodiments of the present invention and are not used to limit the present invention. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principles of the present invention shall be included within the protection scope of the present invention.

Claims

1. An optimized method for suppressing multiple artifacts by second-order blind identification based on an optically pumped magnetometer, characterized in that It includes the following steps: Step 1, collect multi-channel OPM-MEG data, where the multi-channel OPM-MEG data includes brain signal data and artifact data; Step 2, based on the collected data, perform joint approximate diagonalization of the single-time-delay covariance matrix, construct a sensitivity index, evaluate the contribution degree of specific delay points to the separated data, and generate a data-driven parameter adjustment strategy; Step 3, based on the prior knowledge of artifacts, construct artifact templates for low-specificity artifacts and high-specificity artifacts respectively, and form an objective function; Step 4, perform SOBI source separation on the collected multi-channel OPM-MEG data, and optimize the separation process based on the generated parameter adjustment strategy and the objective function to achieve artifact suppression.

2. The second-order blind identification multi-artifact suppression optimization method based on an optically pumped magnetometer according to claim 1, wherein In step 1, a multi-channel OPM sensor array is used to collect brain signal data, and a specific reference sensor is used to collect eye movement, heartbeat, and metal artifact data.

3. The second-order blind identification multi-artifact suppression optimization method based on an optically pumped magnetometer according to claim 2, wherein Denote the data acquisition model as , where represents the data recorded by the OPM sensors, represents the number of OPM sensors, represents the number of sampling points, is the source activity signal, is the number of source activity signals, represents the linear mixing matrix, which describes the mixing mapping from each source activity to the sensor space, represents the random noise.

4. The second-order blind identification multi-artifact suppression optimization method based on an optically pumped magnetometer according to claim 1, wherein Step 2 includes: Step 2.1: Based on the data recorded by the OPM sensor , construct a single-delay covariance matrix , based on the orthogonal transformation matrix transform the single-delay covariance matrix into a diagonal matrix, and the optimal transformation matrix , where τ represents a single time delay, represents the sum of the squares of the non-diagonal elements of the matrix, is the rank transformation; Step 2.2: By constructing sensitivity indicators to quantify the diagonalization trend of a single time-delay matrix and the imbalance of source signal energy distribution to evaluate the contribution of a specific single time delay to the separation of the current data, which is defined as follows: （1） （2） （3） where norm represents the normalization operation, and max(X) represents the maximum element value in matrix X, is the eigenvalue matrix of matrix arranged in descending order of variance. denotes the i-th eigenvalue of the eigenvalue matrix, q is the number of eigenvalues taken in a fixed quantity, the superscript T represents the transpose, m is the total number of eigenvalues, and its value is equal to the number of OPM sensors ; Step 2.

3. For each possible single time delay calculate the sensitivity index , and perform smoothing using Gaussian filtering; adopt a multi-scale peak detection method to identify the sensitivity index peaks at different time scales, and take the delay points within a preset proportion around the peaks to form peak clusters; use the DBSCAN clustering algorithm to identify different peak clusters and group them into intervals, and the set of intervals is denoted as , where K represents the total number of intervals; Step 2.4, based on the divided interval set, construct a data-driven parameter adjustment strategy; considering the sensitivity difference of the source signal to different time delays, introduce a selection probability P for selecting different delay intervals; introduce an interval weight W to assign different weights to intervals containing different time scales; combine the interval selection probability P and the interval weight W, select multiple intervals and perform random sampling of delay points respectively, thereby constituting the optimization constraint conditions for parameter adaptive adjustment and forming a complete parameter adjustment strategy.

5. The second-order blind identification multi-artifact suppression optimization method based on an optically pumped magnetometer according to claim 4, characterized in that, The constraint conditions of step 2.4 are as follows: （4） （5） （6） （7） Among them, represents the number of discrete points included in the j-th interval, q represents the power parameter of the non-linear attenuation function, and represents the set of available intervals when selecting the t-th delay point, represents the selected interval when selecting the t-th delay point, and "\” represents the set difference operation, represents the selected t-th delay point, and U represents the uniform distribution.

6. The second-order blind identification multi-artifact suppression optimization method based on an optically pumped magnetometer according to claim 1, characterized in that, Step 3 includes: Step 3.1, classify artifacts into low-specificity artifact signals with relatively fixed signal characteristics or propagation patterns that can be characterized, and high-specificity artifact signals that require relying on reference sensors due to large individual differences or influencing factors; Step 3.2, for low-specificity artifact signals, establish corresponding time, frequency, and spatial domain prior templates, cluster and compress the scale of the MEG sensor array using the prior knowledge of artifacts, match the decomposition results of local data with the prior templates, and identify the most similar components; if the similarity is above the threshold, update the template for subsequent global identification; if the similarity is below the threshold, it indicates that the artifact characteristics are not significant and are not considered; Step 3.3, for high-specificity artifacts, use their specific quantization indicators and combine the reference channel information, and use mutual information to extract deep correlations for identification and score calculation; Step 3.4, calculate the matching scores for the identified potential various artifact components; （8） （9） （10） （11） Among them, and represent the matching scores of low and high specificity artifacts respectively, represents the Pearson correlation coefficient, where x and y represent the blind separation component source signal and the template signal respectively, represents the fast Fourier transform, represents the absolute value operation, where a and b are the corresponding spatial activation patterns respectively, represents the mutual information index, represents a specific quantization index constructed according to artifact analysis, and represent the i-th values of the blind separation component source signal x and the template signal y respectively, and are the means of the blind separation component source signal x and the template signal y respectively, n is the number of sampling points, p(x,y) represents the joint probability distribution, and p(x) and p(y) are the marginal probability distributions of the blind separation component source signal x and the template signal y respectively, and are specific values of the blind separation component source signal x and the template signal y respectively; Step 3.5, maximizing the objective function of separating known artifacts, that is, the optimization problem is: （12） （13） Among them, and respectively represent the matching scores of the e-th low-specificity artifact and the f-th high-specificity artifact under a set of k delay parameters. is the penalty coefficient; N is the number of finally eliminated components; Nmax is the upper threshold of the eliminated components. represents the total set of optional delays.

7. The second-order blind identification multi-artifact suppression optimization method based on an optically pumped magnetometer according to claim 1, characterized in that In step 4, first whiten the collected OPM-MEG data B to obtain a whitened data matrix Z; calculate the sensitivity index of the data matrix Z at different time delay points, and perform interval division according to the index values to determine the initial delay point set; Based on the set of initial delay points, the SOBI method is used to perform blind separation on the data matrix Z, and the orthogonal transformation matrix V is optimized and solved through Jacobi iteration to obtain the separated source signal components and their corresponding spatial activation patterns; combined with the objective function constructed in step 3, matching recognition and quantitative evaluation of artifact components are performed on the separation results; combined with the parameter adjustment strategy generated in step 2, a directed dynamic reallocation of the delay combinations is performed; a feedback mechanism is established through an intelligent optimization algorithm to adjust and optimize the selection of the time-delay combinations to improve the separation effect and quality in a multi-artifact interference environment.

8. Second-order blind identification multi-artifact suppression optimization system based on an optically pumped magnetometer, characterized in that, Comprising: A data acquisition module for acquiring multi-channel OPM-MEG data, where the multi-channel OPM-MEG data includes brain signal data and artifact data; A strategy construction module for jointly approximately diagonalizing the single-time-delay covariance matrix based on the acquired data, constructing a sensitivity index to evaluate the contribution degree of specific delay points to the separated data, and generating a data-driven parameter adjustment strategy; An objective function construction module for constructing artifact templates for low-specificity artifacts and high-specificity artifacts respectively based on artifact prior knowledge to form an objective function; An optimization module for performing SOBI source separation on the acquired multi-channel OPM-MEG data, and optimizing the separation process based on the generated parameter adjustment strategy and the objective function to achieve artifact suppression.

9. An electronic device, characterized in that, Comprising: One or more processors; A memory for storing one or more programs; Wherein, when the one or more programs are executed by the one or more processors, the one or more processors implement the second-order blind identification multi-artifact suppression optimization method based on an optically pumped magnetometer according to any one of claims 1-7.

10. A computer-readable storage medium, characterized in that, Stored thereon are executable instructions that, when executed by a processor, enable the processor to implement the second-order blind identification multi-artifact suppression optimization method based on an optically pumped magnetometer according to any one of claims 1-7.

Images