CN112461224A - Magnetometer calibration method based on known attitude angle - Google Patents

Magnetometer calibration method based on known attitude angle Download PDF

Info

Publication number
CN112461224A
CN112461224A CN202011247295.4A CN202011247295A CN112461224A CN 112461224 A CN112461224 A CN 112461224A CN 202011247295 A CN202011247295 A CN 202011247295A CN 112461224 A CN112461224 A CN 112461224A
Authority
CN
China
Prior art keywords
magnetometer
coordinate system
magnetic field
axis
matrix
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202011247295.4A
Other languages
Chinese (zh)
Other versions
CN112461224B (en
Inventor
旷俭
李泰宇
牛小骥
刘韬
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Wuhan University WHU
Original Assignee
Wuhan University WHU
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Wuhan University WHU filed Critical Wuhan University WHU
Priority to CN202011247295.4A priority Critical patent/CN112461224B/en
Publication of CN112461224A publication Critical patent/CN112461224A/en
Application granted granted Critical
Publication of CN112461224B publication Critical patent/CN112461224B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C17/00—Compasses; Devices for ascertaining true or magnetic north for navigation or surveying purposes
    • G01C17/38—Testing, calibrating, or compensating of compasses
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C25/00—Manufacturing, calibrating, cleaning, or repairing instruments or devices referred to in the other groups of this subclass
    • G01C25/005—Manufacturing, calibrating, cleaning, or repairing instruments or devices referred to in the other groups of this subclass initial alignment, calibration or starting-up of inertial devices
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R35/00—Testing or calibrating of apparatus covered by the other groups of this subclass

Abstract

The invention provides a magnetometer calibration method based on a known attitude angle. Based on the assumption that the total magnetic field strength of a fixed position is unchanged and the projection component of the magnetic field strength of the fixed position in the same coordinate system is unchanged, the observed value of the magnetometer is projected to a local coordinate system by using a known attitude angle according to a magnetometer measurement model, and magnetometer calibration parameters are obtained through least square iterative solution. Compared with the prior calibration technology, the method has the advantages of simple algorithm and higher efficiency. Meanwhile, the calibration of the soft and hard magnetic effects of the magnetometer is realized, and the coordinate system of the magnetometer is aligned with the coordinate system of the inertial sensor, so that the information of the magnetometer and the inertial sensor is put in a common frame, and the magnetometer and the inertial sensor can play an important role in the fusion positioning of the multisource sensor and the sensor array.

Description

Magnetometer calibration method based on known attitude angle
Technical Field
The invention belongs to the field of magnetometer calibration, and particularly relates to a magnetometer calibration method when an attitude angle is known.
Background
Magnetometers are sensors used to measure magnetic induction, which can provide a reference to magnetic north, a key source of information for attitude estimation in low cost high performance navigation systems. The magnetometer has the problems of zero offset error, cross-axis coupling error, scale factor error and the like due to the limitation of a processing process, and the measurement precision of the magnetometer is influenced; magnetometers are furthermore susceptible to magnetic disturbances of the surroundings, which disturbances are classified as soft magnetic disturbances and hard magnetic disturbances. Therefore, in order to measure more accurate magnetic field information, the magnetometer needs to be calibrated. In addition, in the multi-source fusion sensor positioning, the information of the magnetometer and the inertial sensor is put in a common frame to be the basis for realizing high-precision positioning, but due to the installation process, soft magnetic interference and other reasons, the shafting of the magnetometer and the inertial sensor is difficult to be completely consistent, so that the cross alignment calibration of the magnetometer and the inertial sensor is a key point and the problem to be solved urgently is solved.
At present, according to domestic and foreign documents, ellipsoid fitting based on a maximum likelihood estimation method is the most common magnetometer internal calibration method at present, the calibration of the method is independent of the attitude, and the fact that the measurement value of a magnetometer in a local position or a uniform magnetic field is constant regardless of the direction is utilized. For cross calibration between sensors, this is done mainly by relying on local gravity information, but this calibration method relies on accelerometer measurements.
In summary, when calibrating a magnetometer, how to quickly complete the internal calibration and cross calibration of the magnetometer is a key problem to be solved. The invention provides a simple and effective magnetometer calibration and inertial sensor cross calibration method by using a least square method, including but not limited to the application scene, without depending on attitude and gravity information.
Disclosure of Invention
Aiming at the requirements of rapid calibration of the magnetometer and the problem of misalignment of the magnetometer axis and the IMU inertial sensor axis, the invention provides a method for completing the intrinsic calibration of the magnetometer and the cross calibration of the IMU inertial sensor based on the least square principle under the condition that the attitude angle is known, so that the cross-axis coupling error, the scale factor error, the zero offset error, the soft magnetic error and the hard magnetic error of the calibrated magnetometer can be rapidly realized, and the alignment of the magnetometer and the IMU inertial sensor axis is realized. The method has the advantages of no need of any external equipment or parameter setting, simplicity, feasibility, good universality, capability of meeting the requirement of on-site rapid (about 20 s) calibration and good calibration precision.
The invention adopts the following technical scheme: a method based on least square accomplishes the internal calibration of magnetometer and aligns with the axis system of IMU inertial sensor, finish calibrating mainly through rotating the smartphone, under the presumption that the strength of geomagnetic field is invariable in triaxial magnetic field component and total intensity of magnetic field under the local coordinate system based on the intensity of certain position, through modeling to the measurement model of the magnetometer, use Taylor's expansion and least square method, get the calibration parameter of the magnetometer; the technical scheme comprises the following steps:
step 1, selecting a position in an outdoor open area, and inquiring an earth magnetic field reference model according to longitude and latitude coordinates of the position to obtain a magnetic field total intensity reference value of the position;
step 2, taking the measurement center of the sensor as a rotation center, enabling the sensor to rotate for multiple circles around X, Y, Z three axes of a carrier coordinate system, and acquiring attitude angles (namely pitch angle, roll angle and course angle) and magnetic field information of each epoch in the calibration process;
and 3, establishing a magnetometer measurement model, determining parameters to be estimated, carrying out Taylor expansion on the magnetometer measurement model, and iteratively solving the parameters to be estimated by using a least square method.
Furthermore, when modeling the measurement model of the magnetometer, the measurement model of the magnetometer in the carrier coordinate system is as follows in consideration of the soft magnetic error, the scale factor error, the quadrature axis coupling error, the zero offset error, the hard iron error and the sensor noise of the magnetometer:
Figure BDA0002770454800000021
for magnetometer raw observation data mb_initMean filtering is carried out to reduce the noise n of the magnetometerMTo obtain mb(ii) a Will be provided with
Figure BDA0002770454800000022
The product of (A) is used as a matrix to be estimated and is recorded as
Figure BDA0002770454800000023
B is toMAnd bHIThe sum is taken as a matrix to be estimated and is marked as BMThe model is simplified as follows, and the parameter to be estimated is
Figure BDA0002770454800000024
And BM；
Figure BDA0002770454800000025
The system b is a carrier coordinate system b, the system b takes an inertial sensor IMU phase center as an origin, an x axis points to the advancing direction of the carrier, a y axis is perpendicular to the x axis and points to the right side of the carrier, and a z axis is perpendicular to the x axis and the y axis and forms a right-hand system;
Figure BDA0002770454800000026
representing the actual value of the magnetic field intensity under the carrier coordinate system;
Figure BDA0002770454800000027
an attitude rotation matrix representing the coordinate system of the magnetometer to the coordinate system of the inertial sensor IMU is a 3 multiplied by 3 square matrix; cSIRepresenting a soft magnetic effect change matrix; cNOThe three axes of the magnetometer are converted from non-orthogonality to orthogonality, and the change matrix is a 3 multiplied by 3 square matrix, and the main diagonal element is 0; sMIs a scale factor, is a 3 multiplied by 3 square matrix, only the main diagonal elements have data, and the other elements are 0; m isb_initRepresenting raw observation data of a magnetometer; bMRepresents magnetometer zero offset; bHIRepresenting hard magnetic effect errors; n isMIs the magnetometer noise vector.
Further, when iterative computation is performed on the magnetometer measurement model, first m is usedbAs reference true value of the magnetic field strength in the carrier coordinate system during the first iteration calculation, i.e.
Figure BDA0002770454800000028
Then when the calculation is iterated
Figure BDA0002770454800000029
By passing
Figure BDA00027704548000000210
Then, the true value of the magnetic field intensity in the local coordinate system is obtained
Figure BDA00027704548000000211
The relationship is calculated as follows:
Figure BDA0002770454800000031
in the formula, k represents the kth epoch, and j represents the total number of epochs;
Figure BDA0002770454800000032
an attitude rotation matrix representing the k epoch carrier coordinate system to the local coordinate system,
Figure BDA0002770454800000033
representing a magnetic field intensity reference true value under a k epoch carrier coordinate system; n represents a local coordinate n system, the n system takes the IMU phase center of the inertial sensor as an origin, the x axis is parallel to the local horizontal plane and points to the true north, the y axis is parallel to the local horizontal plane and points to the true east, the z axis is vertical to the local horizontal plane and points downwards, and the three form a right-handed system;
under the assumption that the total strength is not changed based on the magnetic strength, according to Mn_CGRFTo obtain the total magnetic field strength Mn_CGRFBy calculating
Figure BDA0002770454800000034
And Mn_CGRFThe proportional relation of | |, will
Figure BDA0002770454800000035
The data of each axis is changed in equal proportion to obtain Mn_tureTake it as a local seatThe magnetic field truth value under the coordinate system is obtained, and the magnetic field strength reference truth value under the carrier coordinate system is further obtained
Figure BDA0002770454800000036
Mn_tureAnd
Figure BDA0002770454800000037
the calculation process is as follows:
Figure BDA0002770454800000038
Figure BDA0002770454800000039
wherein the content of the first and second substances,
Figure BDA00027704548000000310
representing an attitude rotation matrix from a local coordinate system to a carrier coordinate system;
further, the model will be simplified by a first order Taylor expansion
Figure BDA00027704548000000311
Linearization, namely iteratively solving parameters to be estimated according to a least square method, wherein the specific implementation mode is as follows:
parameter to be estimated
Figure BDA00027704548000000312
Is a 3 x 3 square matrix, and the parameter B to be estimatedMIs a 3 × 1 matrix, and is specifically represented as:
Figure BDA00027704548000000313
general formula
Figure BDA00027704548000000314
The material is spread out and then is put into a bag,
Figure BDA00027704548000000315
the total number of the parameters to be solved is 12, which is marked as X,
X＝[a1 a2 a3 a4 a5 a6 a7 a8 a9 b1 b2 b3]T
because the equation is nonlinear, the equation is developed according to Taylor formula, a second-order term and a high-order term are discarded, and the superscript i represents the iteration number, which is as follows:
Figure BDA0002770454800000041
where Δ represents the residual, let:
Figure BDA0002770454800000042
order: x is the number ofi＝Xi-Xi-1，
Figure BDA0002770454800000043
Equation (7) reduces to:
Figure BDA0002770454800000044
order to
Figure BDA0002770454800000045
Equation (8) is written in the form of a residual equation,
Figure BDA0002770454800000046
the calculation result of the formula (9) can be obtained by the least square method, each observed value is equal in weight, P is a unit array, and the calculation result is obtained
x＝(HTPH)-1HTPz (10)
Xi＝xi+Xi-1 (11)
In the formula, XiFor the raw observation data m of the magnetometerb_inM obtained after mean value filteringbBy passing
Figure BDA0002770454800000051
Updating
Figure BDA0002770454800000052
The iterative process is repeated until convergence, and a parameter result is obtained.
Further, the initial value of the least square method is obtained by the following method,
for the
Figure BDA0002770454800000053
The angle of rotation between the magnetometer and the inertial sensor coordinate system is small, so
Figure BDA0002770454800000054
Approaching a unit array; cSIThe values of all elements are small; cNOThe main diagonal element is 0, and the other elements are smaller; sMOnly main diagonal elements exist, and the numerical value of each element is close to 1; therefore, it is
Figure BDA0002770454800000055
The initial value of (A) is a unit matrix I3×3；
For BMBy further simplifying the model, taking into account only the magnetometer's contribution to zero-bias error, i.e.
Figure BDA0002770454800000056
Raw observation data m of magnetometerb_initAs the true value of the magnetic field in the carrier coordinate system at the first iteration, i.e.
Figure BDA0002770454800000057
According to each epoch
Figure BDA0002770454800000058
(
Figure BDA0002770454800000059
An attitude rotation matrix representing the vector coordinate system of each epoch to the local coordinate system), and the mean value of the magnetic field intensity of each axis in the local coordinate system is obtained
Figure BDA00027704548000000510
As a true value of the magnetic field strength in the local coordinate system. Under the assumption that the total strength is not changed based on the magnetic strength, according to Mn_CGRFTo obtain the total magnetic field strength Mn_CGRFL. By calculating
Figure BDA00027704548000000511
And Mn_CGRFThe proportional relation of | |, will
Figure BDA00027704548000000512
The data of each axis is changed in equal proportion to obtain Mn_tureThis is taken as the true value of the magnetic field in the local coordinate system. Further obtaining the reference true value of the magnetic field intensity under the carrier coordinate system
Figure BDA00027704548000000513
The original data mb_initAnd
Figure BDA00027704548000000514
the difference is calculated to obtain the zero offset of each epoch, and then the zero offset is averaged to obtain BM；
Figure BDA00027704548000000515
When B is obtained in this iterationMAnd B obtained from the last iterationMWhen the difference of each axis is less than 0.1mGauss, the zero offset calibration is considered to be completed, otherwise, the updating is carried out
Figure BDA00027704548000000516
Repeating the above process to obtain BMThis is taken as an initial value in the taylor expansion.
The invention has the following beneficial effects:
(1) the invention completes the internal calibration of the magnetometer and the cross calibration of the inertial sensor on the basis of obtaining the high-precision attitude angle, can reduce zero offset error, cross-axis coupling error, scale factor error, soft magnetic interference and hard magnetic interference and shafting misalignment error of the magnetometer and the inertial sensor, and has better calibration effect.
(2) The method is simple to operate, easy to realize and short in operation time. The calibration work can be completed only by taking the center of the smart phone (namely the carrier) as a rotation center and rotating 720 degrees around X, Y, Z three axes of a carrier coordinate system, the posture is not required, and the calibration efficiency is high.
(3) The invention realizes the alignment of the magnetometer coordinate system and the inertial sensor coordinate system, so that the information of the magnetometer and the inertial sensor is put in a common frame, and the invention can play an important role in the sensor array.
Drawings
Fig. 1 is a schematic diagram of a calibration operation of a mobile phone.
FIG. 2 is a flow chart for iteratively solving magnetometer calibration parameters.
Detailed Description
The technical solution of the present invention is further explained with reference to the drawings and the embodiments.
Step 1, selecting an outdoor open position as a calibration position, acquiring longitude and latitude coordinates and height of the position, and acquiring magnetic field intensity M of the position by inquiring an international Geomagnetic Reference field IGRF (international geographic Reference field) modeln_IGRF；
Step 2, as shown in fig. 1, taking the center of the smart phone (i.e. the carrier) as a rotation center, and rotating the smart phone around X, Y, Z three axes of the carrier coordinate system for 720 degrees;
step 3, obtaining the high precision in the calibration process through pseudo position constraint and reverse smoothing by the method of the patent' MEMS gyroscope automatic calibration methodFrom the carrier coordinate system of each epoch to the attitude rotation matrix of the local coordinate system
Figure BDA0002770454800000061
(the role of the attitude rotation matrix is to effect a projective transformation of the vector from one coordinate system to another);
step 4, establishing a magnetometer measurement model, establishing a relation between the magnetic field strength under a local coordinate system and the magnetic field strength under a carrier coordinate system, obtaining calibration parameters through iterative calculation of the magnetometer measurement model, wherein the calculation flow is shown in fig. 2;
furthermore, the implementation of step 4 comprises the following sub-steps,
41) a measurement model of a magnetometer is modeled. Under the influence of manufacturing process, machining precision and the like, the magnetometer has cross-axis coupling error, scale factor error, zero offset error and sensor noise. In addition, magnetometers are susceptible to magnetic disturbances from the surrounding environment, which are classified as soft and hard magnetic disturbances, which behave in the same manner as quadrature axis coupling errors, scale factor errors, and null offset errors. Therefore, the calibration of the error is required to be completed at the same time. In addition, there is also misalignment error between the magnetometer axis and the inertial sensor axis. The invention considers that a carrier coordinate system is the same as an inertial sensor coordinate system axis system, so that a measurement model of the magnetometer in the carrier coordinate system is as follows:
Figure BDA0002770454800000062
in the formula, b represents a system b of a carrier coordinate system, the system b takes an inertial sensor IMU phase center as an origin, an x axis points to the advancing direction of a carrier, a y axis is perpendicular to the x axis and points to the right side of the carrier, and a z axis is perpendicular to the x axis and the y axis and forms a right-handed system; n represents a local coordinate n system, the n system takes the IMU phase center of the inertial sensor as an origin, the x axis is parallel to the local horizontal plane and points to the true north, the y axis is parallel to the local horizontal plane and points to the true east, the z axis is vertical to the local horizontal plane and points downwards, and the three form a right-handed system;
Figure BDA0002770454800000071
representing the actual value of the magnetic field intensity under the carrier coordinate system;
Figure BDA0002770454800000072
an attitude rotation matrix representing the coordinate system of the magnetometer to the coordinate system of the inertial sensor IMU is a 3 multiplied by 3 square matrix; cSIRepresenting a soft magnetic effect change matrix; cNOThe three axes of the magnetometer are converted from non-orthogonality to orthogonality, and the change matrix is a 3 multiplied by 3 square matrix, and the main diagonal element is 0; sMIs a scale factor, is a 3 multiplied by 3 square matrix, only the main diagonal elements have data, and the other elements are 0; m isb_initRepresenting raw observation data of a magnetometer, bMRepresents magnetometer zero offset; bHIRepresenting hard magnetic effect errors; n isMIs the magnetometer noise vector.
42) For magnetometer raw observation data mb_initMean filtering is carried out to reduce the noise n of the magnetometerMTo obtain mb。
43) In the formula (1)
Figure BDA0002770454800000073
The product of (A) is used as a matrix to be estimated and is recorded as
Figure BDA0002770454800000074
By mbInstead of mb_initAnd nMSumming; b is toMAnd bHIThe sum is taken as a matrix to be estimated and is marked as BM. The magnetometer measurement model is simplified by the following formula, and the parameter to be estimated is
Figure BDA0002770454800000075
And BM。
Figure BDA0002770454800000076
44) In mbAs reference true value of the magnetic field strength in the carrier coordinate system during the first iteration calculation, i.e.
Figure BDA0002770454800000077
Then when the calculation is iterated
Figure BDA0002770454800000078
By passing
Figure BDA0002770454800000079
Then, the true value of the magnetic field intensity in the local coordinate system is obtained
Figure BDA00027704548000000710
The relationship is calculated as follows:
Figure BDA00027704548000000711
wherein k (k ═ 1,2.. j) represents the kth epoch, and j represents the total number of epochs;
Figure BDA00027704548000000712
an attitude rotation matrix representing the k epoch carrier coordinate system to the local coordinate system,
Figure BDA00027704548000000713
and representing the reference true value of the magnetic field intensity in the k epoch carrier coordinate system.
45) Under the assumption that the total strength is not changed based on the magnetic strength, according to Mn_CGRFTo obtain the total magnetic field strength Mn_CGRFBy calculating
Figure BDA00027704548000000714
And Mn_CGRFThe proportional relation of | |, will
Figure BDA00027704548000000715
The data of each axis is changed in equal proportion to obtain Mn_ture. The magnetic field intensity is used as a magnetic field true value under a local coordinate system, and a magnetic field intensity reference true value under a carrier coordinate system is further obtained
Figure BDA00027704548000000716
Mn_tureAnd
Figure BDA00027704548000000717
the calculation process is as follows:
Figure BDA0002770454800000081
Figure BDA0002770454800000082
wherein the content of the first and second substances,
Figure BDA0002770454800000083
representing an attitude rotation matrix from a local coordinate system to a carrier coordinate system;
46) performing first-order Taylor expansion on the simplified model (formula 2), completing the linearization treatment of the model, and then solving the parameter to be estimated by using a least square method
Figure BDA0002770454800000084
And BMThe concrete implementation mode is as follows,
parameter to be estimated
Figure BDA0002770454800000085
Is a 3 x 3 square matrix, and the parameter B to be estimatedMIs a 3 × 1 matrix, and is specifically represented as:
Figure BDA0002770454800000086
general formula
Figure BDA0002770454800000087
The material is spread out and then is put into a bag,
Figure BDA0002770454800000088
the total number of the parameters to be solved is 12, which is marked as X,
X＝[a1 a2 a3 a4 a5 a6 a7 a8 a9 b1 b2 b3]T
because the equation is nonlinear, the equation is developed according to Taylor formula, a second-order term and a high-order term are discarded, and the superscript i represents the iteration number, which is as follows:
Figure BDA0002770454800000089
where Δ represents the residual, let:
Figure BDA0002770454800000091
order: x is the number ofi＝Xi-Xi-1，
Figure BDA0002770454800000092
Equation (7) reduces to:
Figure BDA0002770454800000093
order to
Figure BDA0002770454800000094
Equation (8) is written in the form of a residual equation,
Figure BDA0002770454800000095
the calculation result of the formula (9) can be obtained by a least square method, and the invention considers that all the observed values are equally weighted and P is a unit array to obtain
x＝(HTPH)-1HTPz (10)
Xi＝xi+Xi-1 (11)
XiIs the value of m in step 42)b(ii) a By passing
Figure BDA0002770454800000096
Updating
Figure BDA0002770454800000097
And repeating the steps 43) to 46) until the iteration converges to obtain the parameter result to be estimated.
Further, the initial value of the least square method is obtained by,
due to the fact that
Figure BDA0002770454800000098
The angle of rotation between the magnetometer and the inertial sensor coordinate system is small, so
Figure BDA0002770454800000099
Near unit array, CSIThe values of the respective elements are small, CNOMajor diagonal element is 0, the remaining elements are smaller, SMOnly the main diagonal elements, each element having a value close to 1, are present, so
Figure BDA0002770454800000101
The initial value of (A) is a unit matrix I3×3；
For BMThe initial values are given as follows:
for the mobile phone magnetometer, the error caused by the zero offset is much larger than the errors of the cross-axis coupling, the scale factor and the like, so that when taylor expansion is performed to obtain an initial value, only the influence of the zero offset error on the magnetometer is considered, and the measurement model of the magnetometer is further simplified as follows:
Figure BDA0002770454800000102
raw observation data m of magnetometerb_initAs the true value of the magnetic field in the carrier coordinate system at the first iteration, i.e.
Figure BDA0002770454800000103
According to each epoch
Figure BDA0002770454800000104
The average value of the magnetic field intensity of each axis in the local coordinate system is obtained and used as the true value of the magnetic field intensity in the local coordinate system, the data of j epochs are set, and the calculation process is as follows:
Figure BDA0002770454800000105
under the assumption that the total strength is not changed based on the magnetic strength, according to Mn_CGRFTo obtain the total magnetic field strength Mn_CGRFBy calculating
Figure BDA0002770454800000106
And Mn_CGRFThe proportional relation of | |, will
Figure BDA0002770454800000107
The data of each axis is changed in equal proportion to obtain Mn_tureThe magnetic field intensity is used as a magnetic field true value under the local coordinate system, and a magnetic field intensity reference true value under the carrier coordinate system is further obtained
Figure BDA0002770454800000108
Mn_tureAnd
Figure BDA0002770454800000109
the calculation process is as follows:
Figure BDA00027704548000001010
according to each epoch
Figure BDA00027704548000001011
The reason why the magnetic field strength in the carrier coordinate system of each epoch is obtained and the magnetic field strength in the carrier coordinate system of all epochs is differentiated under the condition that only the zero offset error is considered is because of the existence of the magnetic field strengthDue to zero offset, because the calibration action is axial rotation, the true value of the magnetic field strength under the carrier coordinate system can be obtained by averaging, and the calculation process is as follows:
Figure BDA00027704548000001012
the original data mb_initAnd in the above formula
Figure BDA00027704548000001013
The difference is calculated to obtain the zero offset of each epoch, and then the zero offset is averaged to obtain BM
Figure BDA00027704548000001014
When B is obtained in this iterationMAnd B obtained from the last iterationMWhen the difference of each axis is less than 0.1mGauss, the zero offset calibration is considered to be completed, otherwise, the zero offset calibration is updated according to the formula (12)
Figure BDA0002770454800000111
The above process is repeated. Finally, find BMThis is taken as an initial value in the taylor expansion.
The specific embodiments described herein are merely illustrative of the spirit of the invention. Various modifications or additions may be made to the described embodiments or alternatives may be employed by those skilled in the art without departing from the spirit or ambit of the invention as defined in the appended claims.

Claims (6)

1. A magnetometer calibration method based on a known attitude angle is characterized by comprising the following steps:
step 1, selecting a position in an outdoor open area, and inquiring an earth magnetic field reference model according to longitude and latitude coordinates of the position to obtain a magnetic field total intensity reference value of the position;
step 2, taking the measurement center of the sensor as a rotation center, enabling the sensor to rotate for a plurality of circles around X, Y, Z three axes of a carrier coordinate system respectively, and acquiring attitude angle and magnetic field information of each epoch in the calibration process;
and 3, establishing a magnetometer measurement model, determining parameters to be estimated, carrying out Taylor expansion on the magnetometer measurement model, and iteratively solving the parameters to be estimated by using a least square method.
2. A magnetometer calibration method based on known attitude angles according to claim 1, characterized in that: when modeling a measurement model of the magnetometer, considering that the magnetometer is subjected to soft magnetic errors, scale factor errors, cross-axis coupling errors, zero offset errors, hard iron errors and sensor noise, the measurement model of the magnetometer in a carrier coordinate system is as follows:
Figure FDA0002770454790000011
for magnetometer raw observation data mb_initMean filtering is carried out to reduce the noise n of the magnetometerMTo obtain mb(ii) a Will be provided with
Figure FDA0002770454790000012
The product of (A) is used as a matrix to be estimated and is recorded as
Figure FDA0002770454790000013
B is toMAnd bHIThe sum is taken as a matrix to be estimated and is marked as BM(ii) a After model simplification, the following formula is shown, and the parameter to be estimated is
Figure FDA0002770454790000014
And BM；
Figure FDA0002770454790000015
Wherein b represents a carrier coordinate system b system, and the b system is the IMU phase of the inertial sensorThe center is an origin point, the x axis points to the advancing direction of the carrier, the y axis is perpendicular to the x axis and points to the right side of the carrier, and the z axis is perpendicular to the x axis and the y axis and forms a right-handed system;
Figure FDA0002770454790000016
representing the actual value of the magnetic field intensity under the carrier coordinate system;
Figure FDA0002770454790000017
an attitude rotation matrix representing the coordinate system of the magnetometer to the coordinate system of the inertial sensor IMU is a 3 multiplied by 3 square matrix; cSIRepresenting a soft magnetic effect change matrix; cNOThe three axes of the magnetometer are converted from non-orthogonality to orthogonality, and the change matrix is a 3 multiplied by 3 square matrix, and the main diagonal element is 0; sMIs a scale factor, is a 3 multiplied by 3 square matrix, only the main diagonal elements have data, and the other elements are 0; m isb_initRepresenting raw observation data of a magnetometer; bMRepresents magnetometer zero offset; bHIRepresenting hard magnetic effect errors; n isMIs the magnetometer noise vector.
3. A magnetometer calibration method based on known attitude angles according to claim 2, characterized in that: when iterative calculation is carried out on the magnetometer measurement model, m is used firstlybAs reference true value of the magnetic field strength in the carrier coordinate system during the first iteration calculation, i.e.
Figure FDA00027704547900000112
Then when the calculation is iterated
Figure FDA00027704547900000111
By passing
Figure FDA00027704547900000110
Then, the true value of the magnetic field intensity in the local coordinate system is obtained
Figure FDA0002770454790000021
The relationship is calculated as follows:
Figure FDA0002770454790000022
in the formula, k represents the kth epoch, and j represents the total number of epochs;
Figure FDA00027704547900000215
an attitude rotation matrix representing the k epoch carrier coordinate system to the local coordinate system,
Figure FDA0002770454790000023
representing a magnetic field intensity reference true value under a k epoch carrier coordinate system; n represents a local coordinate n system, the n system takes the IMU phase center of the inertial sensor as an origin, the x axis is parallel to the local horizontal plane and points to the true north, the y axis is parallel to the local horizontal plane and points to the true east, the z axis is vertical to the local horizontal plane and points downwards, and the three form a right-handed system;
under the assumption that the total strength is not changed based on the magnetic strength, according to Mn_CGRFTo obtain the total magnetic field strength Mn_CGRFBy calculating
Figure FDA0002770454790000024
And Mn_CGRFThe proportional relation of | |, will
Figure FDA0002770454790000025
The data of each axis is changed in equal proportion to obtain Mn_tureThe magnetic field intensity is used as a magnetic field true value under the local coordinate system, and a magnetic field intensity reference true value under the carrier coordinate system is further obtained
Figure FDA0002770454790000026
Mn_tureAnd
Figure FDA0002770454790000027
the calculation process is as follows:
Figure FDA0002770454790000028
Figure FDA0002770454790000029
wherein the content of the first and second substances,
Figure FDA00027704547900000210
an attitude rotation matrix representing the local coordinate system to the carrier coordinate system.
4. A magnetometer calibration method based on known attitude angles according to claim 3 characterised in that: simplifying the model by a first order Taylor expansion
Figure FDA00027704547900000211
Linearization, namely iteratively solving parameters to be estimated according to a least square method, wherein the specific implementation mode is as follows:
parameter to be estimated
Figure FDA00027704547900000212
Is a 3 x 3 square matrix, and the parameter B to be estimatedMIs a 3 × 1 matrix, and is specifically represented as:
Figure FDA00027704547900000213
general formula
Figure FDA00027704547900000214
The material is spread out and then is put into a bag,
Figure FDA0002770454790000031
the total number of the parameters to be solved is 12, which is marked as X,
X＝[a1 a2 a3 a4 a5 a6 a7 a8 a9 b1 b2 b3]T
because the equation is nonlinear, the equation is developed according to Taylor formula, a second-order term and a high-order term are discarded, and the superscript i represents the iteration number, which is as follows:
Figure FDA0002770454790000032
where Δ represents the residual, let:
Figure FDA0002770454790000033
order: x is the number ofi＝Xi-Xi-1，
Figure FDA0002770454790000034
Equation (7) reduces to:
Figure FDA0002770454790000035
order to
Figure FDA0002770454790000036
Equation (8) is written in the form of a residual equation,
Figure FDA0002770454790000041
solving the result of equation (9) by least square method, weighting each observed value equally, and obtaining P as unit matrix
x＝(HTPH)-1HTPz (10)
Xi＝xi+Xi-1 (11)
In the formula, XiIs made of magnetismRaw observation data m of force meterb_initM obtained after mean value filteringbBy passing
Figure FDA0002770454790000042
Updating
Figure FDA0002770454790000043
The iterative process is repeated until convergence, and a parameter result is obtained.
5. A magnetometer calibration method based on known attitude angles according to claim 4, characterized in that: the initial value of the least square method is obtained by the following method,
for the
Figure FDA0002770454790000044
The angle of rotation between the magnetometer and the inertial sensor coordinate system is small, so
Figure FDA0002770454790000045
Near unit array, CSIThe values of the respective elements are small, CNOMajor diagonal element is 0, the remaining elements are smaller, SMOnly the main diagonal elements, each element having a value close to 1, are present, so
Figure FDA0002770454790000046
The initial value of (A) is a unit matrix I3×3；
For BMBy further simplifying the model, taking into account only the magnetometer's contribution to zero-bias error, i.e.
Figure FDA0002770454790000047
Raw observation data m of magnetometerb_initAs the true value of the magnetic field in the carrier coordinate system at the first iteration, i.e.
Figure FDA0002770454790000048
According to each epoch
Figure FDA0002770454790000049
Figure FDA00027704547900000410
An attitude rotation matrix from the carrier coordinate system to the local coordinate system for each epoch is expressed, and the average value of the magnetic field intensity of each axis in the local coordinate system is obtained
Figure FDA00027704547900000411
As a true value of the magnetic field strength in the local coordinate system; under the assumption that the total strength is not changed based on the magnetic strength, according to Mn_CGRFTo obtain the total magnetic field strength Mn_CGRFBy calculating
Figure FDA00027704547900000412
And Mn_CGRFThe proportional relation of | |, will
Figure FDA00027704547900000413
The data of each axis is changed in equal proportion to obtain Mn_tureThe magnetic field intensity is used as a magnetic field true value under the local coordinate system, and a magnetic field intensity reference true value under the carrier coordinate system is further obtained
Figure FDA00027704547900000414
The original data mb_initAnd
Figure FDA0002770454790000051
the difference is calculated to obtain the zero offset of each epoch, and then the zero offset is averaged to obtain BM；
Figure FDA0002770454790000052
When B is obtained in this iterationMAnd B obtained from the last iterationMWhen the difference of each axis is less than 0.1mGauss, the zero offset calibration is considered to be completed, otherwise, the zero offset calibration is carried out according to the condition
Figure FDA0002770454790000053
Updating
Figure FDA0002770454790000054
Repeating the above process to obtain BMThis is taken as an initial value in the taylor expansion.
6. A magnetometer calibration method based on known attitude angles according to claim 1, characterized in that: in step 2, the sensor is rotated 720 degrees around X, Y, Z three axes of the carrier coordinate system, and the attitude angle comprises a pitch angle, a roll angle and a heading angle.
CN202011247295.4A 2020-11-10 2020-11-10 Magnetometer calibration method based on known attitude angle Active CN112461224B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011247295.4A CN112461224B (en) 2020-11-10 2020-11-10 Magnetometer calibration method based on known attitude angle

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011247295.4A CN112461224B (en) 2020-11-10 2020-11-10 Magnetometer calibration method based on known attitude angle

Publications (2)

Publication Number Publication Date
CN112461224A true CN112461224A (en) 2021-03-09
CN112461224B CN112461224B (en) 2021-09-14

Family

ID=74825540

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011247295.4A Active CN112461224B (en) 2020-11-10 2020-11-10 Magnetometer calibration method based on known attitude angle

Country Status (1)

Country Link
CN (1) CN112461224B (en)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113074752A (en) * 2021-03-11 2021-07-06 清华大学 Dynamic calibration method and system for vehicle-mounted geomagnetic sensor
CN113310505A (en) * 2021-06-15 2021-08-27 苏州挚途科技有限公司 External parameter calibration method and device of sensor system and electronic equipment
CN114383631A (en) * 2021-12-10 2022-04-22 中国兵器工业集团第二一四研究所苏州研发中心 Real-time calibration method based on least square, Taylor expansion and comprehensive residual combination
CN115839726A (en) * 2023-02-23 2023-03-24 湖南二零八先进科技有限公司 Method, system and medium for jointly calibrating magnetic sensor and angular speed sensor

Citations (18)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5165269A (en) * 1990-10-29 1992-11-24 Iimorrow, Inc. Electronic flux gate compass calibration technique
EP1985233A1 (en) * 2007-04-25 2008-10-29 Commissariat à l'Energie Atomique Method and device for detecting a substantially invariant axis of rotation
CN102735268A (en) * 2012-07-10 2012-10-17 中国人民解放军国防科学技术大学 Strapdown three-shaft magnetometer calibrating method based on posture optimization excitation
CN103175616A (en) * 2011-12-20 2013-06-26 弗卢克公司 Thermal imaging camera with compass calibration
DE102012013516A1 (en) * 2012-07-06 2014-01-09 Giesecke & Devrient Gmbh Calibrating a magnetic sensor
CN103630137A (en) * 2013-12-02 2014-03-12 东南大学 Correction method used for attitude and course angles of navigation system
CN104613983A (en) * 2015-02-03 2015-05-13 中国航天时代电子公司 Whole machine magnetometer calibration method applied to micro unmanned plane
CN105043380A (en) * 2015-06-29 2015-11-11 武汉大学 Indoor navigation method based on a micro electro mechanical system, WiFi (Wireless Fidelity) positioning and magnetic field matching
CN105180968A (en) * 2015-09-02 2015-12-23 北京天航华创科技股份有限公司 IMU/magnetometer installation misalignment angle online filter calibration method
CN105571578A (en) * 2015-12-14 2016-05-11 武汉大学 In-situ rotating modulating north-seeking method utilizing pseudo-observation instead of precise turntable
CN105676302A (en) * 2015-11-12 2016-06-15 东南大学 Magnetometer random noise error compensation method based on improved least square method
CN106323334A (en) * 2015-06-25 2017-01-11 中国科学院上海高等研究院 Magnetometer calibration method based on particle swarm optimization
CN108225378A (en) * 2018-01-25 2018-06-29 陕西土豆数据科技有限公司 A kind of computational methods of compass and accelerometer fix error angle
CN109459020A (en) * 2018-12-24 2019-03-12 哈尔滨工程大学 A kind of inertia and magnetometer combine self-adapting anti-jamming method
US20190133531A1 (en) * 2017-02-13 2019-05-09 National Tsing Hua University Object pose measurement system based on mems imu and method thereof
CN110174123A (en) * 2019-05-08 2019-08-27 苏州大学 A kind of Magnetic Sensor real-time calibration method
EP3620752A1 (en) * 2018-09-05 2020-03-11 Sword Health, S.A. Method and system for calibrating sensors
CN112334736A (en) * 2018-06-13 2021-02-05 西斯纳维 Method for calibrating a magnetometer of an object

Patent Citations (18)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5165269A (en) * 1990-10-29 1992-11-24 Iimorrow, Inc. Electronic flux gate compass calibration technique
EP1985233A1 (en) * 2007-04-25 2008-10-29 Commissariat à l'Energie Atomique Method and device for detecting a substantially invariant axis of rotation
CN103175616A (en) * 2011-12-20 2013-06-26 弗卢克公司 Thermal imaging camera with compass calibration
DE102012013516A1 (en) * 2012-07-06 2014-01-09 Giesecke & Devrient Gmbh Calibrating a magnetic sensor
CN102735268A (en) * 2012-07-10 2012-10-17 中国人民解放军国防科学技术大学 Strapdown three-shaft magnetometer calibrating method based on posture optimization excitation
CN103630137A (en) * 2013-12-02 2014-03-12 东南大学 Correction method used for attitude and course angles of navigation system
CN104613983A (en) * 2015-02-03 2015-05-13 中国航天时代电子公司 Whole machine magnetometer calibration method applied to micro unmanned plane
CN106323334A (en) * 2015-06-25 2017-01-11 中国科学院上海高等研究院 Magnetometer calibration method based on particle swarm optimization
CN105043380A (en) * 2015-06-29 2015-11-11 武汉大学 Indoor navigation method based on a micro electro mechanical system, WiFi (Wireless Fidelity) positioning and magnetic field matching
CN105180968A (en) * 2015-09-02 2015-12-23 北京天航华创科技股份有限公司 IMU/magnetometer installation misalignment angle online filter calibration method
CN105676302A (en) * 2015-11-12 2016-06-15 东南大学 Magnetometer random noise error compensation method based on improved least square method
CN105571578A (en) * 2015-12-14 2016-05-11 武汉大学 In-situ rotating modulating north-seeking method utilizing pseudo-observation instead of precise turntable
US20190133531A1 (en) * 2017-02-13 2019-05-09 National Tsing Hua University Object pose measurement system based on mems imu and method thereof
CN108225378A (en) * 2018-01-25 2018-06-29 陕西土豆数据科技有限公司 A kind of computational methods of compass and accelerometer fix error angle
CN112334736A (en) * 2018-06-13 2021-02-05 西斯纳维 Method for calibrating a magnetometer of an object
EP3620752A1 (en) * 2018-09-05 2020-03-11 Sword Health, S.A. Method and system for calibrating sensors
CN109459020A (en) * 2018-12-24 2019-03-12 哈尔滨工程大学 A kind of inertia and magnetometer combine self-adapting anti-jamming method
CN110174123A (en) * 2019-05-08 2019-08-27 苏州大学 A kind of Magnetic Sensor real-time calibration method

Non-Patent Citations (8)

* Cited by examiner, † Cited by third party
Title
GEBRE-EGZIABHER, D ET AL.: "A gyro-free quaternion-based attitude determination system suitable for implementation using low cast sensors", 《IEEE 2000. POSITION LOCATION AND NAVIGATION SYMPOSIUM》 *
HENKEL, P.ET AL.: "Calibration of magnetic field sensors with two mass-market GNSS receivers", 《2014 11TH WORKSHOP ON POSITIONING, NAVIGATION AND COMMUNICATION (WPNC)》 *
JIAN KUANG ET AL.: "Robust pedestrain dead reckoning based on MEMS-IMU for smartphones", 《SENSORS》 *
JIAN KUANG ET AL.: "Robust Pedestrian Dead Reckoning Based on MEMS-IMU for Smartphones", 《SENSORS》 *
ZHIJIAN DING ET AL.: "Novel low cost calibration methods for MEMS inertial/magnetic integrated sensors", 《PROCEEDINGS OF 2014 IEEE CHINESE GUIDANCE, NAVIGATION AND CONTROL CONFERENCE》 *
张鹏 等: "基于PDR、WiFi指纹识别、磁场匹配组合的室内行人导航定位", 《测绘地理信息》 *
李泰宇 等: "基于磁强计阵列的室内行人定位算法研究", 《传感技术学报》 *
郑元勋 等: "高精度磁信标中心位置与姿态角标定方法", 《中国惯性技术学报》 *

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113074752A (en) * 2021-03-11 2021-07-06 清华大学 Dynamic calibration method and system for vehicle-mounted geomagnetic sensor
CN113310505A (en) * 2021-06-15 2021-08-27 苏州挚途科技有限公司 External parameter calibration method and device of sensor system and electronic equipment
CN113310505B (en) * 2021-06-15 2024-04-09 苏州挚途科技有限公司 External parameter calibration method and device of sensor system and electronic equipment
CN114383631A (en) * 2021-12-10 2022-04-22 中国兵器工业集团第二一四研究所苏州研发中心 Real-time calibration method based on least square, Taylor expansion and comprehensive residual combination
CN115839726A (en) * 2023-02-23 2023-03-24 湖南二零八先进科技有限公司 Method, system and medium for jointly calibrating magnetic sensor and angular speed sensor

Also Published As

Publication number Publication date
CN112461224B (en) 2021-09-14

Similar Documents

Publication Publication Date Title
CN112461224B (en) Magnetometer calibration method based on known attitude angle
EP3411725B1 (en) A method and device for calibration of a three-axis magnetometer
Pei et al. Optimal heading estimation based multidimensional particle filter for pedestrian indoor positioning
CN103153790B (en) The measurement data of the magnetometer using motion sensor and be attached to device estimates equipment and the method for this device yaw angle in gravitational frame of reference
CN106052685B (en) A kind of posture and course estimation method of two-stage separation fusion
CN105785477B (en) The earth magnetism vector measurement error calibrating method that a kind of component combines with total amount constraint
CN105180968A (en) IMU/magnetometer installation misalignment angle online filter calibration method
CN110174123B (en) Real-time calibration method for magnetic sensor
Wahdan et al. Three-dimensional magnetometer calibration with small space coverage for pedestrians
CN106709222B (en) IMU drift compensation method based on monocular vision
CN115507849B (en) Magnetic sensor correction method and system based on INS/GNSS combined navigation assistance
CN104374405A (en) MEMS strapdown inertial navigation initial alignment method based on adaptive central difference Kalman filtering
CN107316280B (en) Li Island satellite image RPC model high-precision geometry location method
CN109000639B (en) Attitude estimation method and device of multiplicative error quaternion geomagnetic tensor field auxiliary gyroscope
CN116817896B (en) Gesture resolving method based on extended Kalman filtering
Liu et al. Tightly coupled modeling and reliable fusion strategy for polarization-based attitude and heading reference system
CN105737850B (en) Mutative scale one direction gravity sample vector matching locating method based on particle filter
CN112577518A (en) Inertial measurement unit calibration method and device
CN112461262A (en) Device and method for correcting errors of three-axis magnetometer
CN106595669B (en) Method for resolving attitude of rotating body
CN113866688B (en) Error calibration method for three-axis magnetic sensor under condition of small attitude angle
CN110398702B (en) Real-time online magnetic calibration method based on multi-sensor fusion
Cui et al. Three-axis magnetometer calibration based on optimal ellipsoidal fitting under constraint condition for pedestrian positioning system using foot-mounted inertial sensor/magnetometer
CN110375773B (en) Attitude initialization method for MEMS inertial navigation system
Chen et al. A novel calibration method for tri-axial magnetometers based on an expanded error model and a two-step total least square algorithm

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant