Magnetic field generator calibration method of medical electromagnetic positioning system

Abstract

The invention discloses a magnetic field generator calibration method of a medical electromagnetic positioning system. The method comprises the following steps: detecting the sensitivity of a single-axis receiving coil; an actual induced electromotive force data set used for calibration is collected, and the induced electromotive force of the single-axis receiving coil under the known pose is collected by building a standard coordinate system frame; a theoretical induced electromotive force data set for calibration is obtained by constructing a magnetic field model; an objective function of iterative optimization is constructed through the two data sets, the placement pose and the equivalent radius of a transmitting coil are optimized by using a Levenberg-Marquardt optimization algorithm, and inherent errors caused by technological processing in a magnetic field generator are calibrated. Under the condition that the induction capacity of the single-axis receiving coil is obtained in advance, independent calibration of the placement poses and the equivalent radiuses of all the transmitting coils in the magnetic field generator is achieved.

Description

Technical Field

[0001] The invention relates to the technical field of medical electromagnetic positioning systems, and in particular to a magnetic field generator calibration method for medical electromagnetic positioning systems. Background Art

[0002] The medical electromagnetic positioning system generally includes a magnetic field generator, a system control unit, a sensor interface unit, several sensors and a computer. A five-degree-of-freedom sensor is a single-axis receiving coil. The determination of the posture of the single-axis receiving coil depends on the difference between the actual induced electromotive force collected and the theoretical induced electromotive force derived from the magnetic field model. This difference is used as the objective function to be optimized. The actual positioning error will be affected by many factors, including the accuracy of the magnetic dipole model, environmental metal interference and the accuracy of the manufacturing process. At present, the calibration methods of medical electromagnetic positioning systems are mostly based on the posture calibration of the receiving sensor. Common methods include global interpolation methods represented by high-order polynomial fitting and local interpolation methods represented by trilinear interpolation and spline interpolation. Nevertheless, the posture-based calibration method is a comprehensive calibration of the accuracy of the entire system, covering all error factors, and there will be problems with improper calibration. The magnetic field model is based on the position and posture of the coil of the magnetic field generator and the inner and outer diameters of the coil itself. It is difficult to ensure the consistency of the size of each excitation coil during the manufacturing process, and the fixture used to fix the coil will also be affected by some manufacturing errors. Therefore, it is necessary to perform a separate calibration of the magnetic field generator before the final posture calibration, and calibrate the equivalent radius and spatial placement of the coil at the same time. Although there are related methods for the calibration of the magnetic field generator, the existing methods do not separate the inductive capability of the single-axis receiving coil from other parameters in the calibration, which will result in the calibration results not being universal when replacing single-axis receiving coils with different inductive capabilities. In addition, the existing magnetic field generator calibration methods do not take into account the processing errors of the coil itself. Therefore, improvements to these methods are expected to improve the overall tracking accuracy of the system, and the calibration results of the transmitting coil can also be used for the positioning of the six-degree-of-freedom receiving coil. Summary of the invention

[0003] The purpose of the present invention is to overcome the shortcomings of the prior art and provide a magnetic field generator calibration method for a medical electromagnetic positioning system. By obtaining the inductive capability of a single-axis receiving coil in advance, the placement posture and equivalent radius of each transmitting coil in the magnetic field generator can be individually calibrated.

[0004] To achieve the above-mentioned purpose, the technical solution provided by the present invention is: a method for calibrating a magnetic field generator of a medical electromagnetic positioning system, wherein the electromagnetic positioning system comprises a magnetic field generator, a system control unit, a sensor interface unit, a plurality of sensors and a computer, wherein one of the five-degree-of-freedom sensors is a single-axis receiving coil; the calibration method only needs to be applied to a single-axis receiving coil, and comprises the following steps:

[0005] S1. Collect the induced electromotive force of the single-axis receiving coil under multiple groups of uniform magnetic fields of different strengths, and fit the relationship between the induced electromotive force and the uniform magnetic field strength, which is called the sensitivity of the single-axis receiving coil; by constructing a magnetic field model, analyze the size of the theoretical magnetic field strength of the single-axis receiving coil under the position when the n transmitting coils in different positions in the magnetic field generator are working, and convert it into theoretical induced electromotive force in combination with sensitivity, analyze the theoretical induced electromotive force of the single-axis receiving coil under different specified positions, and form a theoretical induced electromotive force data set for calibration;

[0006] S2, by time-sharing exciting n transmitting coils in different positions in the magnetic field generator, collecting n induced electromotive forces of the corresponding transmitting coils of the single-axis receiving coil in different specified positions as actual induced electromotive forces, forming an actual induced electromotive force data set for calibration;

[0007] S3. The error vector and objective function of iterative optimization are constructed based on the difference between the corresponding theoretical induced electromotive force and the actual induced electromotive force in the theoretical induced electromotive force data set and the actual induced electromotive force data set. The Levenberg-Marquardt optimization algorithm is used to perform iterative optimization with the posture and equivalent radius of the transmitting coil as the target optimization parameters to obtain the optimal posture and equivalent radius of all transmitting coils, and the offset of the placement posture and equivalent radius is compensated to achieve the calibration of the magnetic field generator.

[0008] Further, in step S1, a tightly wound multi-layer multi-turn coil is processed, and the magnetic field strength B at the center of the multi-layer multi-turn coil is analyzed. c With the excitation current I c Relationship: B c =cI c , c is a parameter describing the linear relationship between the excitation current and the magnetic field strength, wherein the deviation between the simulation and the actual measurement results at the center of the multi-layer multi-turn coil is compared to verify the precision of the processing;

[0009] Taking advantage of the characteristic that a tightly wound multi-layer multi-turn coil generates a uniform magnetic field on the central axis and its vicinity, a single-axis receiving coil is fixed at the center of the multi-layer multi-turn coil to collect the induced electromotive force of the single-axis receiving coil under different field strengths, forming a magnetic field composed of the actual field strength B i and the induced electromotive force E i A dataset composed of i∈[1,Nc ] represents the i-th group of data, N c represents the total number of groups of data;

[0010] When the single-axis receiving coil is fixed at the center of a tightly wound multi-layer multi-turn coil and the magnetic moment directions of the two are consistent, the magnetic field strength B c The relationship between the induced electromotive force E and the corresponding induced electromotive force E is: E = μωSN r B c , where μ is the magnetic induction coefficient, ω is the angular frequency of the excitation current, S and N r are the cross-sectional area and number of turns of the single-axis receiving coil. Except for ω, the other three variables are difficult to measure accurately, so K = μωSN is used. r It represents the inductive capability of the single-axis receiving coil, also known as sensitivity, and K is fitted using the above acquired data set.

[0011] Further, in step S1, a magnetic field model is constructed based on the magnetic dipole model. In the Cartesian coordinate system of the electromagnetic positioning system, the placement coordinates of the jth transmitting coil are w O j =(a j ,b j ,c j ), where a j , b j 、c j They are w O j Coordinate values ​​on the x, y, and z axes, magnetic moment direction vector w H j =(m j ,n j ,p j ), where m j 、n j 、p j They are w H j The components in the x, y, and z axes, and Construct the spherical coordinate system of the jth transmitting coil, with the origin set to w O j , using the coordinates in the spherical coordinate system The Cartesian coordinate system representing the electromagnetic positioning system w H j , where ρ j =1 means w H j The radial distance, θ j and denote the polar angle and azimuth respectively, then and So the basic parameters of the calibration of a single transmitting coil in a magnetic field generator are

[0012] In the Cartesian coordinate system of the electromagnetic positioning system, based on the magnetic dipole model, the position in space is w P=( w x, w y, w z), and the magnetic field simulation is performed by a single-axis receiving coil with an attitude Euler angle of (α, β, γ), where w x、 w y、 w z are w The coordinate values ​​of P on the x, y, and z axes, α, β, and γ are the pitch angle, yaw angle, and roll angle in the Euler angle respectively; the position of the electromagnetic positioning system in the Cartesian coordinate system w The magnetic field vector of P w The calculation formula of B is as follows:

[0013]

[0014] Where, μ0=4π×10 -7 is the magnetic permeability in vacuum, N is the number of turns of the transmitting coil, I g is the amplitude of the exciting current of the transmitting coil, r is the equivalent radius of the magnetic dipole model, and w r= w P- w O j Represents the vector from the center of the transmitting coil to the single-axis receiving coil; Substituting the corresponding parameters, the position of the electromagnetic positioning system in the Cartesian coordinate system can be obtained w The magnetic field magnitude wB of the x, y, and z axes of P x ,wB y ,wB z The calculation method is:

[0015]

[0016] In the formula, the intermediate variable The rotation matrix that converts the Euler angle of the single-axis receiving coil into the Cartesian coordinate system of the electromagnetic positioning system to the Cartesian coordinate system of the single-axis receiving coil The calculation of is as follows:

[0017]

[0018] Thus, the magnetic field of the three axes x, y, and z in the Cartesian coordinate system of the single-axis receiving coil can be obtained. s B x , s B y , s Bz The calculation formula is as follows:

[0019]

[0020] Taking the magnetic moment direction of the uniaxial receiving coil as the z-axis direction of the coordinate system, only s B z Can be converted into the theoretical induced electromotive force ε of the single-axis receiving coil j , with size ε j =K s B z .

[0021] Further, in step S2, a standard coordinate system framework is constructed and a single-axis receiving coil fixture is designed to collect the electromotive force amplitude of the single-axis receiving coil when n different transmitting coils are working, which is called the actual induced electromotive force. In the case of time-sharing excitation of the transmitting coil, the n electromotive force amplitude of the single-axis receiving coil in the mth position is (ε 1m ',ε 2m ',ε 3m ',...,ε jm ',...,ε nm '), where ε nm ' represents the nth electromotive force amplitude in the mth position, and a total of M groups are collected.

[0022] Further, in step S3, when the jth transmitting coil is working, the parameter to be optimized is where r j is the equivalent radius of the transmitting coil as a magnetic dipole model. The theoretical induced electromotive force of the single-axis receiving coil in the mth position is expressed as ε jm (p j ), and the corresponding actual induced electromotive force is ε jm ', select M positions in the tracking domain, then the error vector e between the theoretical induced electromotive force and the actual induced electromotive force j (p j )for:

[0023] e j (p j )=ε j (p j )-ε j '

[0024] In the formula, ε j (p j ) and ε j ' are respectively composed of M ε jm (p j ) and ε jm 'The column vector composed of the objective function to be optimized F cj (pj ) is calculated as follows:

[0025]

[0026] The Levenberg-Marquardt optimization algorithm is used to optimize the parameter p j To optimize, the fixed theoretical position of the transmitting coil and the theoretical inner diameter of the machining are selected as the initial values ​​of the iteration, and p k represents the parameter p in the kth iteration j , the calculation steps of each round of iteration are as follows:

[0027] a. Let the error vector e under the kth iteration j (p k )=(e j1 ,e j2 ,...,e jM ),e jM Yes j (p k ) in the Mth vector element, calculate the Jacobian matrix J of the error vector k and the Heather matrix H k as follows:

[0028]

[0029] In the formula, J k T is the matrix J k The transpose of

[0030] b. Levenberg-Marquardt optimization algorithm combines Gauss-Newton method and steepest descent method, and uses damping coefficient to dynamically adjust the influence of steepest descent method on the iteration step size during the iteration process. The step size of each iteration is d k for:

[0031] d k =(J k T J k +λ k I)- 1 (-J k T ) j (p k )

[0032] Where I is the sixth-order unit matrix, λ k is the damping coefficient;

[0033] c. The determination of the damping coefficient includes the determination of the initial value and the dynamic change during the iteration process. According to the initial Heiser matrix, the initial damping coefficient can be calculated as:

[0034] λk =τ*max(diag(H k ))*||e j (p k )|| 1.5

[0035] Wherein, τ is a pre-set parameter, diag( k ) indicates H k All elements of the diagonal; during the iteration, for each prediction step of the iteration, the scaling factor is calculated The damping coefficient is modified dynamically according to the size of cj (p k ) and F cj (d k +p k ) are the optimization parameters p k and d k +p k The objective function value when is the matrix d k The transpose of k Reflects the accuracy of the prediction step length. k When r is less than 0, the prediction step direction is opposite to the direction of fitness value decrease. The damping coefficient should be increased quickly and the prediction should be abandoned. k ≥0, the forecast should be accepted and the damping coefficient should increase with r k Dynamically adjust the changes in , and improve the consistency between the iteration step direction and the actual fitness value reduction direction;

[0036] d. When the parameters in the iteration process meet any of the following conditions, the iteration ends and the final result is obtained:

[0037] ① Where e1 is the first-order derivative termination condition;

[0038] ②||d k ||<e2(||p k ||+e2), where e2 is the termination condition of the change step distance;

[0039] ③The number of iterations iter>=e3, where e3 is the maximum number of iterations.

[0040] Compared with the prior art, the present invention has the following advantages and beneficial effects:

[0041] 1. The present invention designs a method for measuring the sensitivity of a single-axis receiving coil, which separates the inductive capability of the single-axis receiving coil from other parameters to be calibrated, helps to achieve a separate calibration of the magnetic field generator, and can be applicable to the tracking of multiple single-axis receiving coils with different inductive capabilities.

[0042] 2. The present invention designs a method that can simultaneously calibrate the placement posture of the transmitting coil and the equivalent radius of the transmitting coil when it is approximated as a magnetic dipole model, which is used to calibrate the inherent errors caused by process processing.

[0043] 3. The present invention adopts the Levenberg-Marquardt optimization algorithm to calibrate the magnetic field generator, which can achieve fast and accurate calibration given a suitable initial value. BRIEF DESCRIPTION OF THE DRAWINGS

[0044] Figure 1 is a flow chart of the method of the present invention.

[0045] Figure 2 This is a simulation diagram of the coil magnetic field used for sensitivity detection.

[0046] Figure 3 This is a schematic diagram of the magnetic field model based on the magnetic dipole model.

[0047] Figure 4 It is a diagram of the fixture and standard coordinate system framework for obtaining the theoretical position of the single-axis receiving coil during calibration.

[0048] Figure 5 This is a detailed schematic diagram of a single-axis receiving coil and its fixture.

[0049] Figure 6 It is a working schematic diagram of the transmitting coil and the single-axis receiving coil. DETAILED DESCRIPTION

[0050] The present invention is further described in detail below in conjunction with embodiments and drawings, but the embodiments of the present invention are not limited thereto.

[0051] like Figures 1 to 6 As shown, this embodiment provides a method for calibrating a magnetic field generator of a medical electromagnetic positioning system, wherein the electromagnetic positioning system includes a magnetic field generator, a system control unit, a sensor interface unit, a plurality of sensors, and a computer, wherein one of the five-degree-of-freedom sensors is a single-axis receiving coil. The calibration method only needs to be applied to a single-axis receiving coil. The method, based on the magnetic dipole model and combined with the Levenberg-Marquardt optimization algorithm, realizes the calibration of a single transmitting coil in the magnetic field generator on the premise of obtaining the sensitivity of the single-axis receiving coil. The specific situation is as follows:

[0052] S1. Collect the induced electromotive force of the single-axis receiving coil under multiple groups of uniform magnetic fields of different strengths, and fit the relationship between the induced electromotive force and the uniform magnetic field strength, which is called the sensitivity of the single-axis receiving coil; by constructing a magnetic field model, analyze the size of the theoretical magnetic field strength of the single-axis receiving coil under the position when the n transmitting coils in different positions in the magnetic field generator are working, and convert the sensitivity into the theoretical induced electromotive force, analyze the theoretical induced electromotive force of the single-axis receiving coil under different specified positions, and form a theoretical induced electromotive force data set for calibration; the specific situation is as follows:

[0053] Processing a tightly wound multi-layer multi-turn coil, such as Figure 2 As shown in the figure, the magnetic induction lines at the center of the coil and its vicinity are similar in color and direction, which can be approximated as a uniform magnetic field, and the field strength at the center can be accurately obtained, and the magnetic field strength B at the center of the multi-layer multi-turn coil can be calculated. c With the excitation current I c Relationship: B c =cI c , c is a parameter describing the linear relationship between the excitation current and the magnetic field strength, wherein the deviation between the simulation and the actual measurement results at the center of the multi-layer multi-turn coil is compared to verify the precision of the processing;

[0054] When the single-axis receiving coil is fixed at the center of a tightly wound multi-layer multi-turn coil and the magnetic moment directions of the two are consistent, the magnetic field strength B c The relationship between the induced electromotive force E and the corresponding induced electromotive force E is: E = μωSN r B c , where μ is the magnetic induction coefficient, ω is the angular frequency of the excitation current, S and N r are the cross-sectional area and number of turns of the single-axis receiving coil. Except for ω, the other three variables are difficult to measure accurately, so K = μωSN is used. r It represents the inductive capability of the single-axis receiving coil, also known as sensitivity, and K is fitted using the above acquired data set.

[0055] The single-axis receiving coil is fixed at the center of the processed coil for sensitivity detection, and excited with currents of different magnitudes to collect the actual magnetic field B generated by different excitation currents. i The induced electromotive force E of the single-axis receiving coil under i size, where i∈[1,N c ] represents the i-th group of data, and N represents the total number of data groups. Perform sensitivity straight line fitting on the collected data set: E i =KB i +b, the corresponding sensitivity K and b are calculated as follows:

[0056]

[0057] Where b is the intercept of the fitted line.

[0058] After obtaining the sensitivity of the single-axis receiving coil, the magnetic field model needs to be constructed. The magnetic field model is constructed based on the magnetic dipole model. Figure 3 As shown, in the Cartesian coordinate system of the electromagnetic positioning system, the placement coordinates of the jth transmitting coil are w O j =(a j ,b j ,c j ), where a j , b j 、c j They are w O j Coordinate values ​​on the x, y, and z axes, magnetic moment direction vector w H j =(m j ,n j ,p j ), where m j 、n j 、p j They are w H j The components in the x, y, and z axes, and Construct the spherical coordinate system of the jth transmitting coil, with the origin set to w O j , using the coordinates in the spherical coordinate system Represents the Cartesian coordinate system of the electromagnetic positioning system w H j , where ρ j =1 means w H j The radial distance, θ j and denote the polar angle and azimuth respectively, then and So the basic parameters of the calibration of a single transmitting coil in a magnetic field generator are

[0059] In the Cartesian coordinate system of the electromagnetic positioning system, based on the magnetic dipole model, the position in space is w P=( w x, w y, w z), and the magnetic field simulation is performed by a single-axis receiving coil with an attitude Euler angle of (α, β, γ), where w x、 w y、 w z are wP is the coordinate value of the three axes x, y, and z. α, β, and γ are the pitch angle, yaw angle, and roll angle in the Euler angle respectively. w Magnetic field vector at position P w The calculation formula of B is as follows:

[0060]

[0061] Where, μ0=4π×10 -7 is the magnetic permeability in vacuum, N is the number of turns of the transmitting coil, I g is the amplitude of the coil excitation current, r is the equivalent radius of the magnetic dipole model, and w r= w P- w O j Represents the vector from the center of the transmitting coil to the single-axis receiving coil; Substituting the corresponding parameters, the position of the electromagnetic positioning system in the Cartesian coordinate system can be obtained w The magnitude of the magnetic field in the x, y, and z axes of P w B x , w B y , w B z The calculation method is:

[0062]

[0063] In the formula, the intermediate variable The rotation matrix that converts the Euler angle of the single-axis receiving coil into the Cartesian coordinate system of the electromagnetic positioning system to the Cartesian coordinate system of the single-axis receiving coil The calculation of is as follows:

[0064]

[0065] Thus, the magnetic field of the three axes x, y, and z in the Cartesian coordinate system of the single-axis receiving coil can be obtained. s B x , s B y , s B z The calculation formula is as follows:

[0066]

[0067] Taking the magnetic moment direction of the uniaxial receiving coil as the z-axis direction of the coordinate system, only s B z Can be converted into the theoretical induced electromotive force ε of the single-axis receiving coil j , with size ε j =K s B z .

[0068] S2. By time-sharingly exciting n transmitting coils in different positions in the magnetic field generator, the n induced electromotive forces of the corresponding transmitting coils of the single-axis receiving coil in different specified positions are collected as actual induced electromotive forces to form an actual induced electromotive force data set for calibration. The specific situation is as follows:

[0069] like Figure 4 As shown in the figure, a standard coordinate system frame is built and a single-axis receiving coil fixture is designed to collect the electromotive force amplitude of the single-axis receiving coil when different transmitting coils are working, which is called the actual induced electromotive force. The fixture is fixed in different holes of the standard coordinate system frame to collect the induced electromotive force data of multiple groups of different single-axis receiving coil placement positions. When the single-axis receiving coil is excited in time-sharing mode, the n electromotive force amplitude of the single-axis receiving coil in the mth position is (ε 1m ',ε 2m ',ε 3m ',...,ε jm ',...,ε nm '), where ε nm ' represents the nth electromotive force amplitude in the mth position, that is, the actual induced electromotive force, and a total of M groups are collected. Figure 5 The orientation of the holes machined in the fixture determines the position of the coil in the standard coordinate system. Figure 6 Shown is a working schematic diagram of 8 transmitting coils and a single-axis receiving coil.

[0070] S3. The error vector and objective function of iterative optimization are constructed based on the difference between the theoretical induced electromotive force and the actual induced electromotive force in the theoretical induced electromotive force data set and the actual induced electromotive force data set. The Levenberg-Marquardt optimization algorithm is used to iteratively optimize the posture and equivalent radius of the transmitting coil as the target optimization parameters to obtain the optimal posture and equivalent radius of all transmitting coils, and compensate for the offset of the placement posture and equivalent radius, so as to realize the calibration of the magnetic field generator. The specific situation is as follows:

[0071] When the jth transmitting coil is working, let the parameters to be optimized where r j is the equivalent radius of the transmitting coil as a magnetic dipole model. The theoretical induced electromotive force of the single-axis receiving coil in the mth position is expressed as ε jm (p j ), and the corresponding actual induced electromotive force is ε jm ', select M positions in the tracking domain, then the error vector ej(pj) between the theoretical induced electromotive force and the actual induced electromotive force is:

[0072] e j (p j )=εj (p j )-ε j '

[0073] In the formula, ε j (p j ) and ε j ' are respectively composed of M ε jm (p j ) and ε jm 'The column vector composed of; the objective function to be optimized F cj (p j ) is calculated as follows:

[0074]

[0075] The Levenberg-Marquardt optimization algorithm is used to optimize the parameter p j To optimize, the fixed theoretical position of the transmitting coil and the theoretical inner diameter of the machining are selected as the initial values ​​of the iteration, and p k represents the parameter p in the kth iteration j , the calculation steps of each round of iteration are as follows:

[0076] a. Let the error vector e under the kth iteration j (p k )=(e j1 ,e j2 ,...,e jM ),e jM Yes j (p k ) in the Mth vector element, calculate the Jacobian matrix J of the error vector k and the Heather matrix H k as follows:

[0077]

[0078] H k =J k T J k

[0079] In the formula, J k T is the matrix J k The transpose of

[0080] b. Levenberg-Marquardt optimization algorithm combines Gauss-Newton method and steepest descent method, and uses damping coefficient to dynamically adjust the influence of steepest descent method on the iteration step size during the iteration process. The step size of each iteration is d k for:

[0081] d k =(J kT J k +λ k I) -1 (-J k T ) j (p k )

[0082] Where I is the sixth-order unit matrix, λ k is the damping coefficient;

[0083] c. The determination of the damping coefficient includes the determination of the initial value and the dynamic change during the iteration process. According to the initial Heiser matrix, the initial damping coefficient can be calculated as:

[0084] λ k =τ*max(diag(H k ))*||e j (p k )|| 1.5

[0085] Wherein, τ is a pre-set parameter, diag( k ) indicates H k All elements of the diagonal; during the iteration, for each prediction step of the iteration, the scaling factor is calculated The damping coefficient is modified dynamically according to the size of cj (p k ) and F cj (d k +p k ) are the optimization parameters p k and d k +p k The objective function value when is the matrix d k The transpose of k Reflects the accuracy of the prediction step length. k When r is less than 0, the prediction step direction is opposite to the direction of fitness value decrease. The damping coefficient should be increased quickly and the prediction should be abandoned. k ≥0, the forecast should be accepted and the damping coefficient should increase with r k Dynamically adjust the changes in , and improve the consistency between the iteration step direction and the actual fitness value reduction direction;

[0086] d. When the parameters in the iteration process meet any of the following conditions, the iteration ends and the final result is obtained:

[0087] ① Where e1 is the first-order derivative termination condition;

[0088] ②||dk || <e2(||p k ||+e2), where e2 is the termination condition of the change step distance;

[0089] ③The number of iterations iter>=e3, where e3 is the maximum number of iterations.

[0090] The above embodiments are preferred implementation modes of the present invention, but the implementation modes of the present invention are not limited to the above embodiments. Any other changes, modifications, substitutions, combinations, and simplifications that do not deviate from the spirit and principles of the present invention should be equivalent replacement methods and are included in the protection scope of the present invention.

Claims

1. A method for calibrating a magnetic field generator of a medical electromagnetic positioning system, the electromagnetic positioning system comprising a magnetic field generator, a system control unit, a sensor interface unit, a plurality of sensors and a computer, wherein one of the five-degree-of-freedom sensors is a single-axis receiving coil; the calibration method only needs to be applied to one single-axis receiving coil, characterized in that: The following steps are involved: S1. Collect the induced electromotive force of the single-axis receiving coil under multiple groups of uniform magnetic fields of different strengths, and fit the relationship between the induced electromotive force and the uniform magnetic field strength, which is called the sensitivity of the single-axis receiving coil; by constructing a magnetic field model, analyze the size of the theoretical magnetic field strength of the single-axis receiving coil under the position when the n transmitting coils in different positions in the magnetic field generator are working, and convert it into theoretical induced electromotive force in combination with sensitivity, analyze the theoretical induced electromotive force of the single-axis receiving coil under different specified positions, and form a theoretical induced electromotive force data set for calibration; S2, by time-sharing exciting n transmitting coils in different positions in the magnetic field generator, collecting n induced electromotive forces of the corresponding transmitting coils of the single-axis receiving coil in different specified positions as actual induced electromotive forces, forming an actual induced electromotive force data set for calibration; S3. The error vector and objective function of iterative optimization are constructed based on the difference between the corresponding theoretical induced electromotive force and the actual induced electromotive force in the theoretical induced electromotive force data set and the actual induced electromotive force data set. The Levenberg-Marquardt optimization algorithm is used to perform iterative optimization with the posture and equivalent radius of the transmitting coil as the target optimization parameters to obtain the optimal posture and equivalent radius of all transmitting coils, and the offset of the placement posture and equivalent radius is compensated to achieve the calibration of the magnetic field generator.

2. The method for calibrating a magnetic field generator of a medical electromagnetic positioning system according to claim 1, characterized in that: In step S1, a tightly wound multi-layer multi-turn coil is processed and the magnetic field strength B at the center of the multi-layer multi-turn coil is analyzed. c With the excitation current I c Relationship: B c =cI c , c is a parameter describing the linear relationship between the excitation current and the magnetic field strength, wherein the deviation between the simulation and the actual measurement results at the center of the multi-layer multi-turn coil is compared to verify the precision of the processing; Taking advantage of the characteristic that a tightly wound multi-layer multi-turn coil generates a uniform magnetic field on the central axis and its vicinity, a single-axis receiving coil is fixed at the center of the multi-layer multi-turn coil to collect the induced electromotive force of the single-axis receiving coil under different field strengths, forming a magnetic field composed of the actual field strength B i and the induced electromotive force E i A dataset composed of i∈[1,N c ] represents the i-th group of data, N c represents the total number of groups of data; When the single-axis receiving coil is fixed at the center of a tightly wound multi-layer multi-turn coil and the magnetic moment directions of the two are consistent, the magnetic field strength B c The relationship between the induced electromotive force E and the corresponding induced electromotive force E is: E = μωSN r B c , where μ is the magnetic induction coefficient, ω is the angular frequency of the excitation current, S and N r are the cross-sectional area and number of turns of the single-axis receiving coil. Except for ω, the other three variables are difficult to measure accurately, so K = μωSN is used. r It represents the inductive capability of the single-axis receiving coil, also known as sensitivity, and K is fitted using the above acquired data set.

3. The method for calibrating a magnetic field generator of a medical electromagnetic positioning system according to claim 2, characterized in that: In step S1, a magnetic field model is constructed based on the magnetic dipole model. In the Cartesian coordinate system of the electromagnetic positioning system, the placement coordinates of the jth transmitting coil are w O j =(a j ,b j ,c j ), where a j , b j 、c j They are w O j Coordinate values ​​on the x, y, and z axes, magnetic moment direction vector w H j =(m j ,n j ,p j ), where m j 、n j 、p j They are w H j The components in the x, y, and z axes, and Construct the spherical coordinate system of the jth transmitting coil, with the origin set to w O j , using the coordinates in the spherical coordinate system The Cartesian coordinate system representing the electromagnetic positioning system w H j , where ρ j =1 means w H j The radial distance, θ j and denote the polar angle and azimuth respectively, then and So the basic parameters of the calibration of a single transmitting coil in a magnetic field generator are In the Cartesian coordinate system of the electromagnetic positioning system, based on the magnetic dipole model, the position in space is w P=( w x, w y, w z), and the magnetic field simulation is performed by a single-axis receiving coil with an attitude Euler angle of (α, β, γ), where w x、 w y、 w z are w The coordinate values ​​of P on the x, y, and z axes, α, β, and γ are the pitch angle, yaw angle, and roll angle in the Euler angle respectively; the position of the electromagnetic positioning system in the Cartesian coordinate system w The magnetic field vector of P w The calculation formula of B is as follows: Where, μ0=4π×10 -7 is the magnetic permeability in vacuum, N is the number of turns of the transmitting coil, I g is the amplitude of the exciting current of the transmitting coil, r is the equivalent radius of the magnetic dipole model, and w r= w P- w O j Represents the vector from the center of the transmitting coil to the single-axis receiving coil; Substituting the corresponding parameters, the position of the electromagnetic positioning system in the Cartesian coordinate system can be obtained w The magnitude of the magnetic field in the x, y, and z axes of P w B x , w B y , w B z The calculation method is: In the formula, the intermediate variable The rotation matrix that converts the Euler angle of the single-axis receiving coil into the Cartesian coordinate system of the electromagnetic positioning system to the Cartesian coordinate system of the single-axis receiving coil The calculation of is as follows: Thus, the magnetic field of the three axes x, y, and z in the Cartesian coordinate system of the single-axis receiving coil can be obtained. s B x , s B y , s B z The calculation formula is as follows: Taking the magnetic moment direction of the uniaxial receiving coil as the z-axis direction of the coordinate system, only s B z Can be converted into the theoretical induced electromotive force ε of the single-axis receiving coil j , with size ε j =K s B z .

4. The method for calibrating a magnetic field generator of a medical electromagnetic positioning system according to claim 3, characterized in that: In step S2, a standard coordinate system is constructed and a single-axis receiving coil fixture is designed to collect the electromotive force amplitude of the single-axis receiving coil when n different transmitting coils are working, which is called the actual induced electromotive force. In the case of time-sharing excitation of the transmitting coil, the electromotive force amplitude of the single-axis receiving coil in the mth position is (ε 1m ',ε 2m ',ε 3m ',...,ε jm ',...,ε nm '), where ε nm ' represents the nth electromotive force amplitude in the mth position, and a total of M groups are collected.

5. The method for calibrating a magnetic field generator of a medical electromagnetic positioning system according to claim 4, characterized in that: In step S3, when the jth transmitting coil is working, the parameter to be optimized is where r j is the equivalent radius of the transmitting coil as a magnetic dipole model. The theoretical induced electromotive force of the single-axis receiving coil in the mth position is expressed as ε jm (p j ), and the corresponding actual induced electromotive force is ε jm ', select M positions in the tracking domain, then the error vector e between the theoretical induced electromotive force and the actual induced electromotive force j (p j )for: e j (p j )=e j (p j )-e j ' In the formula, ε j (p j ) and ε j ' are respectively composed of M ε jm (p j ) and ε jm 'The column vector composed of the objective function to be optimized F cj (p j ) is calculated as follows: The Levenberg-Marquardt optimization algorithm is used to optimize the parameter p j To optimize, the fixed theoretical position of the transmitting coil and the theoretical inner diameter of the machining are selected as the initial values ​​of the iteration, and p k represents the parameter p in the kth iteration j , the calculation steps of each round of iteration are as follows: a. Let the error vector e under the kth iteration j (p k )=(e j1 ,e j2 ,...,e jM ),e jM Yes j (p k ) in the Mth vector element, calculate the Jacobian matrix J of the error vector k and the Heather matrix H k as follows: H k =J k T J k In the formula, J k T is the matrix J k The transpose of b. Levenberg-Marquardt optimization algorithm combines Gauss-Newton method and steepest descent method, and uses damping coefficient to dynamically adjust the influence of steepest descent method on the iteration step size during the iteration process. The step size of each iteration is d k for: d k =(J k T J k +λ k I) -1 (-J k T )and j ( p k ) Where I is the sixth-order unit matrix, λ k is the damping coefficient; c. The determination of the damping coefficient includes the determination of the initial value and the dynamic change during the iteration process. According to the initial Heiser matrix, the initial damping coefficient can be calculated as: l k =τ*max(diag(Η k ))*||e j (p k )|| 1.5 Wherein, τ is a pre-set parameter, diag( k ) indicates H k All elements of the diagonal; during the iteration, for each prediction step of the iteration, the scaling factor is calculated The damping coefficient is modified dynamically according to the size of cj (p k ) and F cj (d k +p k ) are the optimization parameters p k and d k +p k The objective function value when is the matrix d k The transpose of k Reflects the accuracy of the prediction step length. k When r is less than 0, the prediction step direction is opposite to the direction of fitness value decrease. The damping coefficient should be increased quickly and the prediction should be abandoned. k ≥0, the forecast should be accepted and the damping coefficient should increase with r k Dynamically adjust the changes in , and improve the consistency between the iteration step direction and the actual fitness value reduction direction; d. When the parameters in the iteration process meet any of the following conditions, the iteration ends and the final result is obtained: ① Where e1 is the first-order derivative termination condition; ②||d k ||<e2(||p k ||+e2), where e2 is the termination condition of the change step distance; ③The number of iterations iter>=e3, where e3 is the maximum number of iterations.

Images