Optimization design method of magnetocardiogram sensor array

Abstract

The invention provides a magnetocardiogram sensor array optimization design method, which comprises the following steps of: performing CT (Computed Tomography) scanning on the upper body of a testee to generate a personalized heart-double lung-trunk model; using a 3D scanner to obtain a testee trunk model containing the magnetocardiogram acquisition board, and registering the testee trunk model with the CT model; based on a volume conduction model and a heart equivalent current dipole source model, establishing a function relationship between a source parameter and a sensor measurement value by adopting a boundary element method; deriving CRB of current dipole parameter estimation errors, and calculating error measurement; and iteratively selecting an optimal sensor position under the limitation of the number of the sensors to obtain an optimal array layout. According to the invention, the source positioning precision of the magnetocardiogram system can be improved.

Description

Technical Field

[0001] The present invention relates to the field of biomedical signal analysis, and more specifically, to a method for optimizing the design of a magnetocardiometric sensor array. Background Art

[0002] Magnetocardiography (MCG) technology is a non-invasive detection method that obtains information on cardiac electrophysiological activity by measuring the weak magnetic field generated by the heart. This technology has important application value in the diagnosis of heart diseases and surgical planning. The MCG system is usually composed of multiple magnetometers to achieve high-precision measurement of the cardiac magnetic field. However, due to cost and practical operation constraints, the number of sensors is often limited, which directly affects the spatial resolution and source localization accuracy of the system. Existing sensor optimization design methods usually rely on specific source estimation algorithms, it is difficult to establish a unified optimization standard, and there are certain limitations in practical applications.

[0003] In the process of implementing the embodiments of the present invention, the inventors found that there are at least the following problems or defects in the prior art: the optimized design of the sensor layout fails to fully utilize the limited sensor resources, resulting in large errors in the positioning of the cardiac magnetic source, and the existing methods lack a universal performance evaluation standard and are difficult to adapt to the needs of different clinical application scenarios. Summary of the invention

[0004] The present invention provides a method for optimizing the design of a magnetocardiometric sensor array, comprising:

[0005] Step 1: Perform CT scan on the upper body of the subject with markers, segment and mesh the obtained CT image, and generate a personalized mesh CT model of the heart, lungs and trunk of the subject;

[0006] Step 2: Use a 3D scanner to obtain a 3D scan model including a high-density magnetic heart acquisition plate and the subject's torso with markers, and align it with the CT model obtained in step 1 to obtain the overall model after alignment;

[0007] Step 3: Based on the volume conduction model of the heart-lungs-trunk obtained in step 2, combined with the equivalent current dipole source model of the heart and the physical model of the magnetometer, the boundary element method based on weighted residual technology is used to solve the forward problem, and the functional relationship between the heart source parameters and the sensor measurement is established. The heart source parameters include position and torque, and then the measurement model of the dipole source is obtained;

[0008] Step 4: Based on the measurement model obtained in step 3 and combined with the noise model, derive the Cramer-Laur bound of the current dipole parameter estimation error;

[0009] Step 5: Under the given sensor quantity limit, the optimal sensor position is selected by iteration to finally obtain the optimal array layout.

[0010] Furthermore, in step 1, the marking point is attached to the side waist of the subject and used as the side waist marking point for CT scanning, and a scanning model is obtained using a 3D scanner. The side waist marking point eliminates the influence of the ups and downs of the abdomen when a person breathes on the model and is not affected by the position of the acquisition board during scanning.

[0011] Furthermore, in step 3, the process of establishing the measurement model includes the following steps:

[0012] The M mesh vertices of the heart are taken as the source space, and each vertex is regarded as a current dipole;

[0013] The heart-lung-torso is modeled as a volume consisting of K uniform and isotropic volume conductor regions, each of which is composed of a closed surface S i separation;

[0014] The surface integral equation of the electric potential is discretized using the boundary element method based on weighted residual technology, and then converted into an algebraic equation for numerical calculation and solution.

[0015] Furthermore, in step 3, for a current dipole with a moment of q at point p, the magnetic field generated at point r is given by the volume integral:

[0016]

[0017] Among them, μ 0 represents the vacuum magnetic permeability, and dV′ represents the volume differential element.

[0018] Furthermore, in step 3, when the current source is located at the heart, the magnetic field b(r) induced outside the body is given by the surface integral:

[0019]

[0020] Where φ(r′) represents the surface S i The electric potential at point r′ on , dS′ represents the directed differential area element at point r′.

[0021] Furthermore, in step 4, the CRB matrix can be obtained by inverting the Fisher information matrix, where the (a, b)th element of the Fisher information matrix is ​​calculated as:

[0022]

[0023] Where A is the linearized vector of the measurement model, C is the noise covariance matrix, and θ a and θ bare the elements in the parameter vector θ.

[0024] Furthermore, in step 4, the error index is calculated, and the volume of the 95% probability concentration ellipsoid is:

[0025]

[0026] where c is the cube root of the eigenvalue of the 3x3 submatrix in the CRB matrix corresponding to the source position variable multiplied by a constant 2.8.

[0027] Furthermore, in step 4, the equivalent uncertainty radius r 95 for:

[0028]

[0029] Where R is the 3x3 submatrix in the CRB matrix corresponding to the source position variable, and der(R) represents the determinant of R.

[0030] Furthermore, in step 5, the median value of the equivalent uncertainty radius of the entire heart source is calculated as follows:

[0031] r 95，median =median(r 95，1 , r 95，2 , ..., r 95，M )

[0032] where r 95，m represents the equivalent uncertainty radius calculated for the mth dipole.

[0033] Furthermore, in step 5, the median of the heart source equivalent uncertainty radius of all remaining possible sensor positions is calculated in each iteration, and the sensor position that can minimize the median is selected, fixed, and then enters the next iteration until the predetermined number of sensors L is reached.

[0034] The above-mentioned embodiments of the present invention have at least the following beneficial effects: the present invention provides a magnetocardiographic sensor array optimization design method, which can improve the source localization accuracy of the magnetocardiographic system. Through the optimization strategy based on the Cramer-Rao bound (CRB), the method can determine the optimal sensor layout for the magnetocardiographic system with a limited number of sensors, thereby minimizing the source localization error. This optimization can not only improve the spatial resolution and sensitivity of the magnetocardiographic system, but also enhance the flexibility and reliability of the sensor array, making it more suitable for clinical applications.

[0035] In addition, this method does not rely on a specific source estimation algorithm, so it has good versatility and adaptability. It can provide the best position layout for any given number of sensors, thereby achieving optimal magnetocardiographic source positioning in different clinical application scenarios. This optimization design independent of specific algorithms can provide a unified performance evaluation standard for magnetocardiography systems, which will help promote the application and development of magnetocardiography technology in cardiac disease diagnosis and surgical planning. BRIEF DESCRIPTION OF THE DRAWINGS

[0036] The above and other objects, features and advantages of the exemplary embodiments of the present invention will become readily understood by reading the detailed description below with reference to the accompanying drawings. In the accompanying drawings, several embodiments of the present invention are shown in an exemplary and non-limiting manner, in which:

[0037] Figure 1 A schematic flow chart of a method for optimizing the design of a magnetocardiometric sensor array provided in one embodiment of the present invention. DETAILED DESCRIPTION

[0038] The principles and spirit of the present invention will be described below with reference to several exemplary embodiments. It should be understood that these embodiments are provided only to enable those skilled in the art to better understand and implement the present invention, and are not intended to limit the scope of the present invention in any way. On the contrary, these embodiments are provided to make the present invention more thorough and complete, and to fully convey the scope of the present invention to those skilled in the art.

[0039] Those skilled in the art know that the embodiments of the present invention can be implemented as a system, device, apparatus, method or computer program product. Therefore, the present invention can be specifically implemented in the following forms, namely: complete hardware, complete software (including firmware, resident software, microcode, etc.), or a combination of hardware and software.

[0040] It should be noted that any number of elements in the drawings is for illustration rather than limitation, and any naming is only for distinction and does not have any limiting meaning.

[0041] Reference below Figure 1 , Figure 1 FIG. 1 is a flow chart of a method for optimizing the design of a magnetic cardiotonic sensor array according to an embodiment of the present invention. Figure 1 As shown, a method 100 for optimizing the design of a magnetocardiographic sensor array includes:

[0042] Step 1: Perform CT scan on the upper body of the subject with markers, segment and mesh the obtained CT image, and generate a personalized mesh CT model of the heart, lungs and trunk of the subject;

[0043] Step 2: Use a 3D scanner to obtain a 3D scan model including a high-density magnetic heart acquisition plate and the subject's torso with markers, and align it with the CT model obtained in step 1 to obtain the overall model after alignment;

[0044] Step 3: Based on the volume conduction model of the heart-lungs-trunk obtained in step 2, combined with the equivalent current dipole source model of the heart and the physical model of the magnetometer, the boundary element method based on weighted residual technology is used to solve the forward problem, and the functional relationship between the heart source parameters and the sensor measurement is established. The heart source parameters include position and torque, and then the measurement model of the dipole source is obtained;

[0045] Step 4: Based on the measurement model obtained in step 3 and combined with the noise model, derive the Cramer-Laur bound of the current dipole parameter estimation error;

[0046] Step 5: Under the given sensor quantity limit, the optimal sensor position is selected by iteration to finally obtain the optimal array layout.

[0047] It should be noted that magnetocardiography is a technology that uses a magnetometer to measure the changes in the weak magnetic field generated by the heart, which is of great significance for the diagnosis and research of heart disease. In this method, a CT scan of the upper body of the subject with markers is first performed, and the obtained CT image is segmented and gridded to generate a personalized grid CT model of the subject's heart-lungs-trunk. This process is to obtain accurate three-dimensional structural information of the subject's heart and trunk, and provide basic data for subsequent sensor layout.

[0048] Specifically, CT scanning is a technology that uses X-rays to scan the human body from different angles and generates detailed images of the internal structure of the human body through computer processing.

[0049] Furthermore, in the present method, the segmentation and meshing of the CT image refers to separating the heart, lungs and trunk regions in the CT image and converting them into a three-dimensional mesh model. This step usually involves medical image processing techniques, such as image segmentation algorithms and three-dimensional reconstruction techniques. The setting of the marker points is to provide reference points in the CT scan and subsequent 3D scans to ensure accurate registration of the model. These marker points are usually set on the side waist of the subject to reduce the impact of respiratory movement on the model.

[0050] More specifically, when performing a CT scan, a high-resolution scanning mode can be selected to obtain clearer image details. During the image segmentation process, a threshold-based method or a more advanced machine learning algorithm can be used to distinguish different tissue types. Meshing involves converting the segmented image into a three-dimensional model composed of multiple small grid cells, the size and shape of which can be adjusted as needed to balance the accuracy and computational efficiency of the model.

[0051] Furthermore, professional medical image processing software can be used to assist in completing these steps, such as 3DSlicer or Mimics. During the 3D scanning process, optical scanners or laser scanners can be used, which can quickly obtain the surface geometry information of the subject's torso and align it with the CT model to ensure the accuracy of the final model. During the alignment process, methods such as feature matching or iterative closest point (ICP) algorithm can be used to achieve this.

[0052] In some embodiments, in step 1, the marking points are attached to the side waist of the subject and used as side waist marking points for CT scanning, and a scanned model is obtained using a 3D scanner. The side waist marking points eliminate the influence of the ups and downs of the abdomen when a person breathes on the model and are not affected by the position of the acquisition board during scanning.

[0053] It should be noted that this embodiment involves attaching marking points to the side waist of the subject when performing a CT scan. The role of the side waist marking points is to provide a reference when performing a CT scan and subsequent 3D scans to ensure the accuracy and consistency of the model. The setting of the marking points helps to eliminate the impact of the ups and downs of the abdomen when a person breathes on the model, and is not affected by the position of the high-density magnetic heart acquisition board during the scanning process. The high-density magnetic heart acquisition board is a device for placing a magnetometer, which includes a square slot that can fix the sensor and is used to design the optimal layout of the sensor.

[0054] Specifically, the lateral waist markers are usually made of materials that are easy to identify in CT images, such as metal balls or special marking tapes. The size and shape of these markers need to be small enough so as not to cause discomfort to the subject or affect the accuracy of the scan results. When setting the markers, it is necessary to ensure that they are evenly distributed on the lateral waist of the subject to provide sufficient reference information in the subsequent model alignment process. The design of the high-density magnetic heart acquisition board needs to take into account the size and layout requirements of the magnetometer to ensure that each magnetometer can be fixed in the best position. In addition, the square slots on the acquisition board should be designed to be deep enough so that the magnetometer can be placed firmly without affecting the comfort of the subject.

[0055] More specifically, when performing a CT scan, a high-resolution scanning mode can be selected to obtain clearer image details. During the image segmentation process, a threshold-based method or a more advanced machine learning algorithm can be used to distinguish different tissue types. Meshing involves converting the segmented image into a three-dimensional model composed of multiple small grid cells, the size and shape of which can be adjusted as needed to balance the accuracy and computational efficiency of the model.

[0056] Furthermore, professional medical image processing software can be used to assist in completing these steps, such as 3DSlicer or Mimics. During the 3D scanning process, optical scanners or laser scanners can be used, which can quickly obtain the surface geometry information of the subject's torso and align it with the CT model to ensure the accuracy of the final model. During the alignment process, methods such as feature matching or iterative closest point (ICP) algorithm can be used to achieve it.

[0057] In some embodiments, in step 3, the process of establishing the measurement model includes the following steps:

[0058] The M mesh vertices of the heart are taken as the source space, and each vertex is regarded as a current dipole;

[0059] The heart-lung-torso is modeled as a volume consisting of K uniform and isotropic volume conductor regions, each of which is composed of a closed surface S i separation;

[0060] The surface integral equation of the electric potential is discretized using the boundary element method based on weighted residual technology, and then converted into an algebraic equation for numerical calculation and solution.

[0061] It should be noted that the present embodiment involves a volume conduction model based on the heart-lungs-trunk, combined with an equivalent current dipole source model of the heart and a physical model of a magnetometer, and uses a boundary element method based on weighted residual technology to solve the forward problem. The forward problem refers to calculating the distribution of the magnetic field generated by the heart on the body surface when the heart source parameters (such as position and torque) are known. The equivalent current dipole source model is a simplified model of cardiac electrophysiological activity, which equates the electrophysiological activity of the heart to one or more current dipoles for ease of calculation and analysis. The physical model of the magnetometer describes the physical process of the magnetometer measuring the magnetic field.

[0062] Specifically, the volume conduction model of the heart-lungs-trunk is established by considering the heart, lungs, and trunk as volumes composed of multiple uniform and isotropic volume conductor regions. These regions are separated by closed surfaces, and the conductivity of each region can be set according to physiological data. Each vertex in the equivalent current dipole source model is regarded as a current dipole, and its current density can be expressed by a specific mathematical expression. The physical model of the magnetometer involves parameters such as the sensitivity, directionality, and measurement noise of the magnetometer. When using the boundary element method based on weighted residual technology, the surface integral equation of the electric potential needs to be discretized and converted into an algebraic equation for numerical calculation and solution. This process involves the segmentation of closed surfaces, the positioning of nodes, and the approximate representation of the electric potential.

[0063] More specifically, when building a measurement model, you can choose the appropriate mesh density and element type to improve the accuracy and efficiency of the calculation. For example, the heart area can use a finer mesh to capture complex electrophysiological activities, while the trunk area can use a coarser mesh. When performing boundary element method calculations, you can choose the appropriate weight function and residual form to optimize the calculation results.

[0064] Furthermore, the model parameters can be verified and adjusted through experiments to ensure the accuracy and reliability of the model. When dealing with the surface integral equation of the magnetic field, efficient numerical integration methods such as Gaussian integration or Simpson integration can be used to improve the accuracy of the calculation. For the physical model of the magnetometer, its performance changes under different environmental conditions, such as the influence of factors such as temperature and humidity on the measurement results, can be considered, and corresponding compensation can be made in the model.

[0065] In some embodiments, in step 3, for a current dipole with moment q located at point p, the magnetic field generated at point r is given by the volume integral:

[0066]

[0067] Among them, μ 0 represents the vacuum magnetic permeability, and dV′ represents the volume differential element.

[0068] It should be noted that this embodiment involves deriving the Cramér-Rao Bound (CRB) of the current dipole parameter estimation error based on the measurement model in combination with the noise model. CRB is a concept in statistical theory, which is used to describe the minimum variance limit of parameter estimation under certain conditions. In this method, CRB is used to evaluate the accuracy of cardiac source parameter estimation, thereby providing a theoretical basis for the optimization of the sensor array. The noise model is used to describe the statistical characteristics of random errors in the measurement process, including the distribution, mean and variance of noise.

[0069] Specifically, the process of deriving CRB involves the linearization of the measurement model and the calculation of the noise covariance matrix. The linearization of the measurement model is achieved by separating the part of the model that only depends on the source parameters from the part that only involves the geometry and conductivity based on the boundary element method, thereby achieving the linearization of the source parameters. The noise covariance matrix describes the statistical characteristics of the measurement noise. It is usually assumed that the noise is Gaussian white noise, and its covariance matrix is ​​a scalar multiple of the unit matrix.

[0070] Furthermore, when calculating CRB, it is necessary to use the Fisher information matrix, which is a matrix related to the measurement model and the noise model, and its inverse matrix gives the CRB. The calculation of the Fisher information matrix involves the derivation of the Jacobian matrix of the measurement model and the processing of the noise covariance matrix. In addition, the calculation of CRB can also be implemented through numerical methods such as Monte Carlo simulation to handle complex measurement models and noise models.

[0071] More specifically, when deriving the CRB, an appropriate numerical method can be selected to calculate the derivative part of the Fisher information matrix, such as the finite difference method or the automatic differentiation method, to improve the accuracy and efficiency of the calculation. For the noise model, the appropriate noise distribution and covariance matrix can be selected according to the actual measurement environment and equipment characteristics. For example, if the measurement noise has obvious non-Gaussian characteristics, other types of distributions can be selected, such as the Laplace distribution or the t distribution.

[0072] Furthermore, the parameters of the noise model can be estimated through experimental data, such as estimating the noise covariance matrix by measuring the sample covariance matrix of the noise. After calculating the CRB, the CRB can be used to evaluate the impact of different sensor layouts on the source localization accuracy, thereby providing guidance for the optimization of the sensor array. For example, the sensor layout with the smallest CRB can be selected as the optimal layout, or the performance of different layouts can be evaluated by comparing their CRBs. In addition, other performance indicators, such as mean square error (MSE) or bias, can also be considered to assist in the optimization of the sensor array.

[0073] In some embodiments, in step 3, when the current source is located in the heart, the magnetic field b′(r) induced outside the body is given by the surface integral:

[0074]

[0075] Where φ(r′) represents the surface S i The electric potential at point r′ on , dS′ represents the directed differential area element at point r′.

[0076] It should be noted that this embodiment involves the process of iteratively selecting the optimal sensor position under a given sensor number limit to ultimately obtain the optimal array layout. The core of this process is how to iteratively select the best sensor placement position through an algorithm under a limited number of sensors to achieve the optimal magnetocardiographic signal acquisition effect. The iterative selection process is an optimization algorithm that gradually adjusts the sensor position, and each selection can improve the overall performance index (such as CRB) until the preset number of sensors is reached or the performance index is no longer significantly improved.

[0077] Specifically, in the process of implementing the iterative selection of the optimal sensor position, it is first necessary to define a performance evaluation index, such as the determinant of the 3x3 submatrix corresponding to the source position variable in the CRB matrix. This index reflects the impact of the sensor layout on the accuracy of cardiac source localization. Then, from all possible sensor positions, a position is selected that minimizes the performance evaluation index. This process may involve calculating the CRB for each potential position and selecting the best one.

[0078] Furthermore, in each iteration, the remaining possible sensor positions need to be updated and the above selection process is repeated until the preset number of sensors is reached. In this process, the selection of sensor positions is based on the calculation of the median of the equivalent uncertainty radius of the cardiac source, which reflects the central trend of the uncertainty of all cardiac source position estimates under the current sensor layout.

[0079] More specifically, in the iterative selection process, heuristic algorithms or optimization algorithms, such as genetic algorithms, simulated annealing algorithms, or particle swarm optimization algorithms, can be used to search for the optimal sensor position. These algorithms can effectively find solutions with better performance in a larger search space. In addition, the mutual influence between sensors, such as the distance and direction between sensors, can also be considered to avoid mutual interference of signals. In each iteration, in addition to calculating the performance evaluation index, the uniformity of the sensor layout can also be considered to ensure uniform coverage of the entire magnetic cardio signal acquisition area.

[0080] Furthermore, if the number or type of sensors is limited in practical applications, alternative solutions can be considered, such as increasing the sensitivity of the sensors or using different types of sensors to make up for the lack of quantity. In addition, the sensor layout can be verified and adjusted through experiments to ensure that the final sensor array can meet the needs of practical applications.

[0081] In some embodiments, in step 4, the CRB matrix can be obtained by inverting the Fisher information matrix, where the (a, b)th element of the Fisher information matrix is ​​calculated as:

[0082]

[0083] Where A is the linearized vector of the measurement model, C is the noise covariance matrix, and θ a and θ are the elements in the parameter vector θ.

[0084] It should be noted that the present embodiment involves calculating error indicators, in particular the volume of the 95% probability concentration ellipsoid. This volume is a statistical concept used to describe the size of the region where the unbiased estimated mean is located at a probability of 95%. In this context, it is used to quantify the uncertainty of the cardiac source position estimate. Specifically, this volume is calculated based on the Cramer-Laurent bound (CRB) matrix, which provides a lower bound on the parameter estimation error and is an important tool for evaluating the accuracy of the estimation.

[0085] Specifically, when calculating the volume of the 95% probability concentration ellipsoid, we first need to determine the 3x3 submatrix in the CRB matrix that corresponds to the source location variable. The cube root of the eigenvalue of this submatrix multiplied by a constant (usually 2.8) can be used to calculate the radius of the ellipsoid. This constant is based on the characteristics of the normal distribution, which ensures that with a probability of 95%, the estimated mean will fall within this ellipsoid. After calculating the radius of the ellipsoid, the formula for the volume of the ellipsoid (4 / 3πr 3 ) to calculate the volume. The size of this volume directly reflects the uncertainty of the estimated cardiac source position. The smaller the volume, the higher the accuracy of the estimate.

[0086] More specifically, during implementation, numerical methods can be used to calculate the CRB matrix, especially when the model is complex or has many parameters. For example, Monte Carlo simulation can be used to estimate the CRB matrix, and a large number of simulation experiments can be used to approximate the true CRB value. In addition, optimization algorithms can be considered to minimize this volume to find the optimal sensor layout. In actual operation, it may be necessary to compare different sensor layout schemes to select the one that provides the minimum 95% probability concentration ellipsoid volume.

[0087] Furthermore, if computational resources are limited in practical applications, one can consider using an approximate CRB calculation method or adopt a heuristic algorithm to approximate the optimal solution. In addition, the performance of the selected sensor layout can be verified experimentally to ensure that it meets the accuracy requirements in practical applications.

[0088] In some embodiments, in step 4, the error index is calculated, and the volume of the 95% probability concentration ellipsoid is:

[0089]

[0090] where c is the cube root of the eigenvalue of the 3x3 submatrix in the CRB matrix corresponding to the source position variable multiplied by a constant 2.8.

[0091] It should be noted that the present embodiment involves calculating the volume of the 95% probability concentration ellipsoid, which is a statistical measure used to quantify the uncertainty of the cardiac source location estimate. In this context, the 95% probability concentration ellipsoid is a geometric body that contains the 95% probability region of the cardiac source location estimate. The volume of this ellipsoid is determined by the CRB matrix, which provides a lower bound on the parameter estimation error and is an important tool for evaluating the estimation accuracy. Specifically, the calculation of this volume involves the eigenvalues ​​of the 3x3 submatrix in the CRB matrix corresponding to the source location variable.

[0092] Specifically, when calculating the volume of the 95% probability concentration ellipsoid, we first need to extract the 3x3 submatrix associated with the cardiac source location from the CRB matrix. The eigenvalues ​​of this submatrix represent the error estimates in different directions, and their cube roots multiplied by a constant (usually 2.8) can be used to calculate the radius of the ellipsoid. This constant is based on the characteristics of the normal distribution, which ensures that with a probability of 95%, the estimated mean will fall within this ellipsoid. Then, using the formula for the volume of the ellipsoid (4 / 3πr 3 ) to calculate the volume, where r is the radius of the ellipsoid. The size of this volume directly reflects the uncertainty of the estimated cardiac source location. The smaller the volume, the higher the accuracy of the estimate. In practical applications, this volume can be used to compare the performance of different sensor layout schemes and select the one that provides the smallest volume.

[0093] More specifically, during implementation, numerical methods can be used to calculate the CRB matrix, especially when the model is complex or has many parameters. For example, Monte Carlo simulation can be used to estimate the CRB matrix, and a large number of simulation experiments can be used to approximate the true CRB value. In addition, optimization algorithms can be considered to minimize this volume to find the optimal sensor layout. In actual operation, it may be necessary to compare different sensor layout schemes to select the one that provides the minimum 95% probability concentration ellipsoid volume.

[0094] Furthermore, if computational resource constraints are encountered in practical applications, an approximate CRB calculation method can be considered, or a heuristic algorithm can be used to approximate the optimal solution. In addition, the performance of the selected sensor layout can be verified experimentally to ensure that it meets the accuracy requirements in practical applications. During the experimental verification process, the performance of the sensor layout can be evaluated using actual magnetic heart data, or its performance can be predicted using simulated data. Through these methods, it can be ensured that the selected sensor layout can provide high-precision cardiac source location estimation.

[0095] In some embodiments, in step 4, the equivalent uncertainty radius r 95 for:

[0096]

[0097] Where R is the 3x3 submatrix in the CRB matrix corresponding to the source position variable, and det(R) represents the determinant of R.

[0098] It should be noted that the present embodiment involves calculating an equivalent uncertainty radius, which is a metric used to quantify the uncertainty of the cardiac source location estimate. The equivalent uncertainty radius is calculated based on the CRB matrix, which provides a simplified metric for representing the uncertainty of the cardiac source location estimate. Specifically, the equivalent uncertainty radius is calculated by the determinant of the 3x3 submatrix corresponding to the source location variable in the CRB matrix, and this radius corresponds to the radius of a sphere whose volume is equal to the volume of the 95% probability concentration ellipsoid calculated previously.

[0099] Specifically, when calculating the equivalent uncertainty radius, it is first necessary to extract the 3x3 submatrix related to the cardiac source position from the CRB matrix. Then, the determinant of this submatrix is ​​calculated, and the cube root of its square root multiplied by a constant (usually 2.8) can be used to calculate the equivalent uncertainty radius. This constant is based on the characteristics of the normal distribution, which ensures that with a probability of 95%, the estimated mean will fall within a sphere with this radius as the radius. This equivalent uncertainty radius provides an intuitive metric for comparing the performance of different sensor layout schemes and selecting the scheme that provides the smallest radius. In practical applications, this radius can be used to evaluate the accuracy of sensor layout and to optimize it when designing sensor arrays.

[0100] More specifically, during implementation, numerical methods can be used to calculate the CRB matrix, especially when the model is complex or has many parameters. For example, Monte Carlo simulation can be used to estimate the CRB matrix, and a large number of simulation experiments can be used to approximate the true CRB value. In addition, optimization algorithms can be considered to minimize this equivalent uncertainty radius to find the optimal sensor layout.

[0101] Furthermore, in actual operation, it may be necessary to compare different sensor layout schemes to select the one that provides the minimum equivalent uncertainty radius. If computing resource constraints are encountered in actual applications, it is possible to consider using an approximate CRB calculation method or a heuristic algorithm to approximate the optimal solution.

[0102] Furthermore, the performance of the selected sensor layout can be verified experimentally to ensure that it meets the accuracy requirements of practical applications. In the experimental verification process, the performance of the sensor layout can be evaluated using actual magnetic heart data, or its performance can be predicted using simulated data. Through these methods, it can be ensured that the selected sensor layout can provide high-precision cardiac source position estimation.

[0103] In some embodiments, in step 5, the median value of the equivalent uncertainty radius of the entire heart source is calculated as follows:

[0104] r 95，median =median(r 95，1 , r 95，2 ,…,r 95，M )

[0105] where r 95，m represents the equivalent uncertainty radius calculated for the mth dipole.

[0106] It should be noted that the present embodiment involves calculating the median of the equivalent uncertainty radius of the entire cardiac source. The median is determined by analyzing the uncertainty of multiple cardiac source position estimates, providing a comprehensive metric for evaluating the overall performance of the sensor array layout. Specifically, the median is obtained by calculating the equivalent uncertainty radius of each cardiac source position estimate and then taking the median of these radii. This method can reduce the impact of outliers and more accurately reflect the central trend of the sensor array layout, thereby better representing the typical uncertainty level of the cardiac source position estimate, thereby providing a more reliable performance indicator for the optimization of the sensor array.

[0107] Specifically, when calculating the median of the equivalent uncertainty radius of the entire cardiac source, it is first necessary to determine the equivalent uncertainty radius of each cardiac source position estimate. This usually involves analyzing the CRB matrix of each cardiac source position, extracting the 3x3 submatrix corresponding to the source position variable, and calculating the cube root of its determinant. Then, the equivalent uncertainty radii of all cardiac source positions are collected and the median of these radii is calculated. This median can be used as an important indicator to evaluate the performance of the sensor array layout.

[0108] Furthermore, in practical applications, this median value can help researchers or engineers compare different sensor layout schemes and select the scheme that can provide the minimum median value, thereby optimizing the layout of the sensor array.

[0109] More specifically, in the implementation process, numerical methods can be used to calculate the CRB matrix for each cardiac source location, especially when the model is complex or has many parameters. For example, Monte Carlo simulation can be used to estimate the CRB matrix, and a large number of simulation experiments can be used to approximate the actual CRB value. In addition, optimization algorithms can be considered to minimize this median value to find the optimal sensor layout.

[0110] Furthermore, in actual operation, it may be necessary to compare different sensor layout schemes to select the one that provides the minimum median value. If computing resources are limited in actual applications, it is possible to consider using an approximate CRB calculation method or a heuristic algorithm to approximate the optimal solution. In addition, the performance of the selected sensor layout can be verified experimentally to ensure that it meets the accuracy requirements in actual applications.

[0111] Preferably, in the experimental verification process, actual magnetic heart data can be used to evaluate the performance of the sensor layout, or simulated data can be used to predict its performance. Through these methods, it can be ensured that the selected sensor layout can provide high-precision cardiac source position estimation.

[0112] In some embodiments, in step 5, the median of the cardiac source equivalent uncertainty radius of all remaining possible sensor positions is calculated in each iteration, and the sensor position that can minimize the median is selected, fixed, and then enters the next iteration until the predetermined number of sensors L is reached.

[0113] It should be noted that this embodiment involves calculating the median of the equivalent uncertainty radius of the cardiac source of all remaining possible sensor positions in each iteration, and selecting the sensor position that can minimize the median. This process is the core of sensor array optimization, and the purpose is to minimize the overall positioning error by iteratively selecting the best position of the sensor. In this process, the median of the equivalent uncertainty radius is a key performance indicator, which comprehensively reflects the uncertainty level of all cardiac source position estimates under the current sensor layout.

[0114] Specifically, in the process of iteratively selecting sensor locations, we first need to define a performance evaluation metric, namely the median of the equivalent uncertainty radius of the heart source. This metric is obtained by calculating the equivalent uncertainty radius corresponding to each potential sensor location and taking the median of these radii. The equivalent uncertainty radius is calculated based on the CRB matrix, which provides a lower bound on the parameter estimation error.

[0115] Furthermore, in each iteration, the algorithm evaluates all remaining possible sensor positions and selects the sensor position that minimizes the median equivalent uncertainty radius. This process involves a lot of calculations, especially when the sensor position selection space is large. Specific parameter settings, such as the number of iterations, the search range of sensor positions, and the step size, need to be adjusted according to actual conditions to ensure the effectiveness and efficiency of the optimization process.

[0116] More specifically, in the process of iteratively selecting sensor locations, a variety of optimization algorithms, such as genetic algorithms, simulated annealing algorithms, or particle swarm optimization algorithms, can be used to search for optimal sensor locations. These algorithms can effectively find solutions with better performance in a larger search space.

[0117] Preferably, the mutual influence between sensors, such as the distance and direction between sensors, can also be considered to avoid mutual interference of signals. In each iteration, in addition to calculating the performance evaluation index, the uniformity of the sensor layout can also be considered to ensure uniform coverage of the entire magnetic cardio signal acquisition area. If the number or type of sensors is limited in actual applications, alternative solutions can be considered, such as increasing the sensitivity of the sensors or using different types of sensors to make up for the lack of quantity.

[0118] Furthermore, the sensor layout can be verified and adjusted experimentally to ensure that the final sensor array can meet the needs of practical applications. During the experimental verification process, the performance of the sensor layout can be evaluated using actual magnetic heart data, or its performance can be predicted using simulated data. Through these methods, it can be ensured that the selected sensor layout can provide high-precision cardiac source position estimation.

[0119] The above-mentioned embodiments of the present invention have the following beneficial effects: the magnetocardiographic sensor array optimization design method described in the present invention can provide a personalized heart-lung-trunk model for the magnetocardiographic system, and combine the high-density magnetocardiographic acquisition board and 3D scanning technology to achieve accurate model registration, thereby laying a foundation for the optimization design of the sensor array. This method can significantly improve the spatial resolution and sensitivity of the magnetocardiographic system, thereby improving the accuracy of source positioning.

[0120] By implementing the implementation steps of the optimization design method, it is possible to ensure that the optimal sensor position is iteratively selected under a limited number of sensors, and finally the optimal array layout is obtained. This method can not only minimize the error of magnetocardiographic source localization, but also improve the flexibility and reliability of the sensor array, which is especially suitable for clinical applications with a limited number of sensors, and improves the source localization accuracy of the magnetocardiographic system.

[0121] Furthermore, the storage medium of the embodiment of the present application stores program instructions that can implement all the above methods, wherein the program instructions can be stored in the above storage medium in the form of a software product, including several instructions for enabling a computer device (which can be a personal computer, server, or network device, etc.) or a processor to execute all or part of the steps of the methods described in each embodiment of the present application. The aforementioned storage medium includes: various media that can store program codes, such as a USB flash drive, a mobile hard disk, a read-only memory (ROM), a random access memory (RAM), a magnetic disk or an optical disk, or terminal devices such as a computer, a server, a mobile phone, and a tablet.

[0122] The above descriptions are only some preferred embodiments of the present invention and an explanation of the technical principles used. Those skilled in the art should understand that the scope of the invention involved in the embodiments of the present invention is not limited to the technical solutions formed by a specific combination of the above technical features, but should also cover other technical solutions formed by any combination of the above technical features or their equivalent features without departing from the above inventive concept. For example, the above features are replaced with (but not limited to) technical features with similar functions disclosed in the embodiments of the present invention to form a technical solution.

Claims

1. A method for optimizing the design of a magnetocardiographic sensor array, characterized in that: The following steps are involved: Step 1: Perform CT scan on the upper body of the subject with markers, segment and mesh the obtained CT image, and generate a personalized mesh CT model of the heart, lungs and trunk of the subject; Step 2: Use a 3D scanner to obtain a 3D scan model including a high-density magnetic heart acquisition plate and the subject's torso with markers, and align it with the CT model obtained in step 1 to obtain the overall model after alignment; Step 3: Based on the volume conduction model of the heart-lungs-trunk obtained in step 2, combined with the equivalent current dipole source model of the heart and the physical model of the magnetometer, the boundary element method based on weighted residual technology is used to solve the forward problem, and the functional relationship between the heart source parameters and the sensor measurement is established. The heart source parameters include position and torque, and then the measurement model of the dipole source is obtained; Step 4: Based on the measurement model obtained in step 3 and combined with the noise model, derive the Cramer-Laur bound of the current dipole parameter estimation error; Step 5: Under the given sensor quantity limit, the optimal sensor position is selected by iteration to finally obtain the optimal array layout.

2. The method for optimizing the design of a magnetocardiometric sensor array according to claim 1, characterized in that: In step 1, the marking point is attached to the side waist of the subject and used as the side waist marking point for CT scanning, and a scanning model is obtained using a 3D scanner. The side waist marking point eliminates the influence of the ups and downs of the abdomen when a person breathes on the model and is not affected by the position of the acquisition board during scanning.

3. The method for optimizing the design of a magnetocardiometric sensor array according to claim 1, characterized in that: In step 3, the process of establishing the measurement model includes the following steps: The M mesh vertices of the heart are taken as the source space, and each vertex is regarded as a current dipole; The heart-lung-torso is modeled as a volume consisting of K uniform and isotropic volume conductor regions, each of which is composed of a closed surface S i Separation; The surface integral equation of the electric potential is discretized using the boundary element method based on weighted residual technology, and then converted into an algebraic equation for numerical calculation and solution.

4. The method for optimizing the design of a magnetocardiometric sensor array according to claim 3, characterized in that: In step 3, for a current dipole with moment q at point p, the magnetic field generated at point r is given by the volume integral: Where μ0 represents the vacuum magnetic permeability and dV′ represents the volume differential element.

5. The method for optimizing the design of a magnetocardiometric sensor array according to claim 3, characterized in that: In step 3, when the current source is located in the heart, the magnetic field b′(r) induced outside the body is given by the surface integral: Where φ(r′) represents the surface S i The electric potential at point r′ on , dS′ represents the directed differential area element at point r′.

6. The method for optimizing the design of a magnetocardiometric sensor array according to claim 1, characterized in that: In step 4, the CRB matrix can be obtained by inverting the Fisher information matrix, where the (a, b)th element of the Fisher information matrix is ​​calculated as: Where A is the linearized vector of the measurement model, C is the noise covariance matrix, and θ a and θ b are the elements in the parameter vector θ.

7. The method for optimizing the design of a magnetocardiometric sensor array according to claim 1, characterized in that: In step 4, the error index is calculated, and the volume of the 95% probability concentration ellipsoid is: where c is the cube root of the eigenvalue of the 3x3 submatrix in the CRB matrix corresponding to the source position variable multiplied by a constant 2.

8.

8. The method for optimizing the design of a magnetocardiometric sensor array according to claim 1, characterized in that: In step 4, the equivalent uncertainty radius r 95 for: Where R is the 3x3 submatrix in the CRB matrix corresponding to the source position variable, and det(R) represents the determinant of R.

9. The method for optimizing the design of a magnetocardiometric sensor array according to claim 1, characterized in that: In step 5, the median value of the equivalent uncertainty radius of the entire heart source is calculated as follows: r 95,median =median(r 95,1 ,r 95,2 ,…,r 95,M ) where r 95,m represents the equivalent uncertainty radius calculated for the mth dipole.

10. The method for optimizing the design of a magnetocardiometric sensor array according to claim 1, characterized in that: In step 5, the median of the heart source equivalent uncertainty radius of all remaining possible sensor positions is calculated in each iteration, and the sensor position that can minimize the median is selected and fixed before entering the next iteration until the predetermined number of sensors L is reached.

Images