CN106197405A

CN106197405A - Inertia earth magnetism matching locating method under the influence of geomagnetic diurnal change

Info

Publication number: CN106197405A
Application number: CN201610617495.1A
Authority: CN
Inventors: 解伟男; 黄黎平; 李清华; 奚伯齐
Original assignee: Harbin Institute of Technology Shenzhen
Current assignee: Harbin Institute of Technology Shenzhen
Priority date: 2016-08-01
Filing date: 2016-08-01
Publication date: 2016-12-07
Anticipated expiration: 2036-08-01
Also published as: CN106197405B

Abstract

The invention provides the inertia earth magnetism matching locating method under the influence of a kind of geomagnetic diurnal change, read the position measurements a of the point to be matched in N number of moment from inertial navigation system_iAnd b_i, obtain geomagnetic field intensity information I from gaussmeter_i；The position of the point N number of to be matched according to inertial navigation system instruction, reads the reference value I (a of corresponding geomagnetic field intensity from the geomagnetic database prestored respectively_i,b_i), reference value I of the gradient of geomagnetic field intensity_x,iAnd I_y,i；Introduce and initialize longitude and latitude error, course error and geomagnetic diurnal change error；Calculate the increment of longitude and latitude error, the increment of course error and the increment δ M of geomagnetic diurnal change error；Update longitude and latitude error, course error and geomagnetic diurnal change error M；Judge whether to meet and terminate iterated conditional, calculate parameter K and δ K according to the M after updating；Obtain longitude and latitude error, course error and geomagnetic diurnal change error according to iterative computation, acquired results substitution coupling track is i.e. obtained with the relation equation of reference locus and mates track.

Description

Inertial geomagnetic matching positioning method under the influence of geomagnetic diurnal variation

技术领域technical field

本发明涉及一种地磁日变影响下的惯性地磁匹配定位方法，属于惯性地磁匹配定位方法技术领域。The invention relates to an inertial geomagnetic matching positioning method under the influence of the daily change of the geomagnetic field, and belongs to the technical field of inertial geomagnetic matching positioning methods.

背景技术Background technique

导航技术随着现代科技的发展日新月异，总的来看主要有惯性导航、卫星导航、地形匹配导航、天文导航、地磁导航等导航技术。在各式各样的导航技术中，地磁导航以其高度自主性、隐蔽性和无积累误差等优势成为当今导航领域的一大研究热门。惯性地磁匹配是指，载体上安装有惯导系统和磁传感器，载体在飞行过程中，惯导系统输出载体的运动轨迹信息，将此输出轨迹称为“参考轨迹”，由于惯性器件存在漂移，因此参考轨迹与载体运动所经历的真实轨迹存在一定误差；在惯导系统输出载体轨迹的同时，由地磁传感器测量得到载体所经历位置的地磁场信息，然后根据预存于计算机内的地球磁场信息、地磁场实时测量信息以及载体参考轨迹，采用匹配方法，得到载体的匹配轨迹。地磁场包含稳定磁场和变化磁场两部分。前者包括主磁场和异常场，是地磁场的主体，也是地磁匹配所要用到的信息；后者起源于固体地球外的空间电流体系，相当于地磁测量时的干扰。地磁日变属于后者，包括太阳静日变化和太阴日变化，其平均变化幅度约为几纳特到几十纳特。目前，地磁日变并未形成较为统一的数学模型，不同地点的地磁日变曲线有较大的差异，同一地点的地磁日变曲线不同季节也不同、白天和晚上也不同，因此在惯性地磁匹配定位过程中难以进行有效的补偿或剔除。所以，地磁日变能够对惯性地磁匹配定位的精度及可靠性产生严重的影响。With the development of modern science and technology, navigation technology is changing with each passing day. In general, there are mainly navigation technologies such as inertial navigation, satellite navigation, terrain matching navigation, astronomical navigation, and geomagnetic navigation. Among all kinds of navigation technologies, geomagnetic navigation has become a hot topic in the field of navigation due to its advantages of high autonomy, concealment and no accumulation of errors. Inertial geomagnetic matching means that an inertial navigation system and a magnetic sensor are installed on the carrier. During the flight of the carrier, the inertial navigation system outputs the movement trajectory information of the carrier. This output trajectory is called the "reference trajectory". Due to the drift of the inertial device, Therefore, there is a certain error between the reference trajectory and the real trajectory experienced by the carrier movement; while the inertial navigation system outputs the carrier trajectory, the geomagnetic sensor measures the geomagnetic field information of the carrier’s experienced position, and then according to the geomagnetic field information pre-stored in the computer, The real-time measurement information of the geomagnetic field and the reference trajectory of the carrier are used to obtain the matching trajectory of the carrier by using a matching method. The earth's magnetic field includes two parts: a stable magnetic field and a changing magnetic field. The former includes the main magnetic field and anomalous field, which is the main body of the geomagnetic field and is also the information used for geomagnetic matching; the latter originates from the space current system outside the solid earth, which is equivalent to the interference in geomagnetic measurement. The geomagnetic diurnal variation belongs to the latter, including solar quiet diurnal variation and lunar diurnal variation, and its average variation range is about a few nanots to tens of nats. At present, the geomagnetic diurnal variation has not formed a relatively unified mathematical model. The geomagnetic diurnal variation curves of different locations are quite different. The geomagnetic diurnal variation curves of the same location are also different in different seasons, day and night. It is difficult to effectively compensate or eliminate during the positioning process. Therefore, the diurnal variation of the geomagnetic field can have a serious impact on the accuracy and reliability of inertial geomagnetic matching positioning.

基于上述原因，对地磁日变影响下的惯性地磁匹配定位方法的需要就十分迫切。然而由于日变的不确定性、复杂性，对日变十分难于处理，现有地磁导航方法中对这方面的处理非常少。在目前的地磁导航研究中有的文献是利用复杂方法拟合地磁日变曲线、有的文献虽然强调地磁日变影响，但是没能够给出处理地磁日变的方法，直接忽略地磁日变的影响。经文献检索，授权公告号为CN103115624B、授权公告日为2014年12月10日的中国发明专利，公开了“一种基于地磁匹配的地磁日变修正方法”，该专利基于FMI方法对日变进行了拟合，而且该方法需要导航当天以及前后共三天的地磁数据为基础进行计算，不易实际应用。谢仕民等在文章“地磁匹配导航关键技术浅析”中指出了地磁变化场(也即地磁日变)对导航精度产生较大误差，但并没能够提出校正地磁日变的方法。李豫泽等在文章“基于ICCP算法的地磁匹配定位方法”(见《科学计算及信息处理》)中则没有提及地磁日变，忽略了地磁日变的影响，无法保证匹配定位精度。Based on the above reasons, there is an urgent need for an inertial geomagnetic matching positioning method under the influence of geomagnetic diurnal variation. However, due to the uncertainty and complexity of diurnal variation, it is very difficult to deal with diurnal variation, and the existing geomagnetic navigation methods deal with this aspect very little. In the current geomagnetic navigation research, some literatures use complex methods to fit the geomagnetic diurnal variation curve, and some literatures emphasize the influence of geomagnetic diurnal variation, but fail to give a method to deal with geomagnetic diurnal variation, directly ignoring the influence of geomagnetic diurnal variation . After literature search, the Chinese invention patent with the authorized announcement number CN103115624B and the authorized announcement date on December 10, 2014 discloses "a method for correcting the daily variation of geomagnetism based on geomagnetic matching". In addition, this method needs to be calculated based on the geomagnetic data of the day of navigation and three days before and after, which is not easy for practical application. Xie Shimin et al. pointed out in the article "Analysis of the Key Technology of Geomagnetic Matching Navigation" that the geomagnetic variation field (that is, the diurnal variation of geomagnetism) has a large error on the navigation accuracy, but they have not been able to propose a method to correct the diurnal variation of geomagnetism. In the article "Geomagnetic Matching and Positioning Method Based on ICCP Algorithm" (see "Scientific Computing and Information Processing"), Li Yuze et al. did not mention the diurnal change of geomagnetism, ignoring the influence of diurnal change of geomagnetism, and cannot guarantee the matching positioning accuracy.

发明内容Contents of the invention

本发明的目的是为了解决上述现有技术存在的问题，进而提供一种地磁日变影响下的惯性地磁匹配定位方法。The object of the present invention is to solve the above-mentioned problems in the prior art, and further provide an inertial geomagnetic matching positioning method under the influence of geomagnetic diurnal variation.

本发明的目的是通过以下技术方案实现的：The purpose of the present invention is achieved through the following technical solutions:

一种地磁日变影响下的惯性地磁匹配定位方法，步骤如下：An inertial geomagnetic matching positioning method under the influence of geomagnetic diurnal variation, the steps are as follows:

步骤一、从惯导系统读取当前时刻以及前N-1个时刻的待匹配点的位置测量值a_i和b_i，其中a_i表示经度，b_i表示纬度，下标i表示不同时刻，i＝1…N，N为整数且N＞2，i为1表示当前时刻，通过地磁传感器获得当前时刻以及前N-1个时刻的地磁场强度的测量值I_i；Step 1. Read the position measurement values a _i and b _i of the points to be matched at the current moment and the previous N-1 moments from the inertial navigation system, where a _i represents the longitude, b _i represents the latitude, and the subscript i represents different moments, i=1...N, N is an integer and N>2, i is 1 to represent the current moment, and the measured value I _i of the geomagnetic field strength at the current moment and the previous N-1 moments is obtained through the geomagnetic sensor;

步骤二、根据惯导系统指示的N个待匹配点的位置，从预先存储的地磁数据库中分别读取相应的地磁场强度的参考值I(a_i,b_i)、地磁场强度的梯度的参考值I_x,i和I_y,i，其中I_x,i表示地磁场强度在经度方向的梯度在第i点位置上的取值，I_y,i表示地磁场强度在纬度方向的梯度在第i点位置上的取值；Step 2. According to the positions of the N points to be matched indicated by the inertial navigation system, read the reference value I(a _i , b _i ) of the corresponding geomagnetic field strength and the gradient of the geomagnetic field strength respectively from the pre-stored geomagnetic database. Reference values I _x,i and I _y,i , where I _x,i represents the value of the gradient of the geomagnetic field strength in the longitude direction at the i-th point, and I _y,i represents the gradient of the geomagnetic field strength in the latitude direction at The value at the position of the i-th point;

步骤三、引入并初始化经纬度误差、航向误差以及地磁日变误差：Step 3. Introduce and initialize latitude and longitude error, heading error and geomagnetic diurnal variation error:

M＝[Δx Δy α δ]^T＝[0 0 0 0]^T M＝[Δx Δy α δ] ^T ＝[0 0 0 0] ^T

步骤四、根据公式(1)、公式(2)和公式(3)计算迭代参数G、F和H：Step 4, calculate iteration parameters G, F and H according to formula (1), formula (2) and formula (3):

G＝g(M) (1)G=g(M) (1)

F＝f(M) (2)F＝f(M) (2)

H＝F^-1 (3)H＝F ^-1 (3)

其中：in:

g(M)＝[g₁(M) g₂(M) g₃(M) g₄(M)]^T；g(M)=[g ₁ (M) g ₂ (M) g ₃ (M) g ₄ (M)] ^T ;

$f f ((M m)) = = [\begin{matrix} {f f}_{1111} ((M m)) & {f f}_{1212} ((M m)) & {f f}_{1313} ((M m)) & {f f}_{1414} ((M m)) \\ {f f}_{1212} ((M m)) & {f f}_{22 twenty two} ((M m)) & {f f}_{23 twenty three} ((M m)) & {f f}_{24 twenty four} ((M m)) \\ {f f}_{1313} ((M m)) & {f f}_{23 twenty three} ((M m)) & {f f}_{3333} ((M m)) & {f f}_{3434} ((M m)) \\ {f f}_{1414} ((M m)) & {f f}_{24 twenty four} ((M m)) & {f f}_{3434} ((M m)) & {f f}_{4444} ((M m)) \end{matrix}];;$

$\begin{matrix} {g g}_{11} ((M m)) = = {Σ Σ}_{i i = = 11}^{N N} {I I}_{x x,, i i} [[{I I}_{x x,, i i} (({a a}_{i i}^{' '} cos cos α α + + {b b}_{i i}^{' '} sin sin α α - - {a a}_{i i}^{' '} + + Δ Δ x x)) \\ + + {I I}_{y the y,, i i} ((- - {a a}_{i i}^{' '} sin sin α α + + {b b}_{i i}^{' '} cos cos α α - - {b b}_{i i}^{' '} + + Δ Δ y the y)) + + {I I}_{t t,, i i} + + δ δ]] \end{matrix};;$

$\begin{matrix} {g g}_{22} ((M m)) = = {Σ Σ}_{i i = = 11}^{N N} {I I}_{y the y,, i i} [[{I I}_{x x,, i i} (({a a}_{i i}^{' '} cos cos α α + + {b b}_{i i}^{' '} sin sin α α - - {a a}_{i i}^{' '} + + Δ Δ x x)) \\ + + {I I}_{y the y,, i i} ((- - {a a}_{i i}^{' '} sin sin α α + + {b b}_{i i}^{' '} cos cos α α - - {b b}_{i i}^{' '} + + Δ Δ y the y)) + + {I I}_{t t,, i i} + + δ δ]] \end{matrix};;$

$\begin{matrix} {g g}_{33} ((M m)) = = {Σ Σ}_{i i = = 11}^{N N} [[{I I}_{x x,, i i} (({a a}_{i i}^{' '} cos cos α α + + {b b}_{i i}^{' '} sin sin α α - - {a a}_{i i}^{' '} + + Δ Δ x x)) + + {I I}_{y the y,, i i} ((- - {a a}_{i i}^{' '} sin sin α α \\ + + {b b}_{i i}^{' '} cos cos α α - - {b b}_{i i}^{' '} + + Δ Δ y the y)) + + {I I}_{t t,, i i} + + δ δ]] \times \times [[{I I}_{x x,, i i} ((- - {a a}_{i i}^{' '} sin sin α α + + {b b}_{i i}^{' '} cos cos α α)) \\ + + {I I}_{y the y,, i i} ((- - {a a}_{i i}^{' '} cos cos α α - - {b b}_{i i}^{' '} sin sin α α))]] \end{matrix};;$

$\begin{matrix} {g g}_{44} ((M m)) = = {Σ Σ}_{i i = = 11}^{N N} [[{I I}_{x x,, i i} (({a a}_{i i}^{' '} cos cos α α + + {b b}_{i i}^{' '} sin sin α α - - {a a}_{i i}^{' '} + + Δ Δ x x)) \\ + + {I I}_{y the y,, i i} ((- - {a a}_{i i}^{' '} sin sin α α + + {b b}_{i i}^{' '} cos cos α α - - {b b}_{i i}^{' '} + + Δ Δ y the y)) + + {I I}_{t t,, i i} + + δ δ]] \end{matrix};;$

${f f}_{1111} ((M m)) = = {Σ Σ}_{i i = = 11}^{N N} {I I}_{x x,, i i}^{22};;$

${f f}_{1212} ((M m)) = = {Σ Σ}_{i i = = 11}^{N N} {I I}_{x x,, i i} {I I}_{y the y,, i i};;$

${f f}_{1313} ((M m)) = = {Σ Σ}_{i i = = 11}^{N N} {I I}_{x x,, i i} [[{I I}_{x x,, i i} ((- - {a a}_{i i}^{' '} s the s i i n no α α + + {b b}_{i i}^{' '} c c o o s the s α α)) + + {I I}_{y the y,, i i} ((- - {a a}_{i i}^{' '} c c o o s the s α α - - {b b}_{i i}^{' '} s the s i i n no α α))]];;$

${f f}_{1414} ((M m)) = = {Σ Σ}_{i i = = 11}^{N N} {I I}_{x x,, i i};;$

${f f}_{22 twenty two} ((M m)) = = {Σ Σ}_{i i = = 11}^{N N} {I I}_{y the y,, i i}^{22};;$

${f f}_{23 twenty three} ((M m)) = = {Σ Σ}_{i i = = 11}^{N N} {I I}_{y the y,, i i} [[{I I}_{x x,, i i} ((- - {a a}_{i i}^{' '} s the s i i n no α α + + {b b}_{i i}^{' '} c c o o s the s α α)) + + {I I}_{y the y,, i i} ((- - {a a}_{i i}^{' '} c c o o s the s α α - - {b b}_{i i}^{' '} s the s i i n no α α))]];;$

${f f}_{24 twenty four} ((M m)) = = {Σ Σ}_{i i = = 11}^{N N} {I I}_{y the y,, i i};;$

$\begin{matrix} {f f}_{3333} ((M m)) = = {Σ Σ}_{i i = = 11}^{N N} {[[{I I}_{x x,, i i} ((- - {a a}_{i i}^{' '} sin sin α α + + {b b}_{i i}^{' '} cos cos α α)) + + {I I}_{y the y,, i i} ((- - {a a}_{i i}^{' '} cos cos α α - - {b b}_{i i}^{' '} sin sin α α))]]}^{22} \\ + + [[{I I}_{x x,, i i} (({a a}_{i i}^{' '} cos cos α α + + {b b}_{i i}^{' '} sin sin α α - - {a a}_{i i}^{' '} + + Δ Δ x x)) + + {I I}_{y the y,, i i} ((- - {a a}_{i i}^{' '} sin sin α α + + {b b}_{i i}^{' '} cos cos α α \\ - - {b b}_{i i}^{' '} + + Δ Δ y the y)) + + {I I}_{t t,, i i} + + δ δ]] \times \times [[{I I}_{x x,, i i} ((- - {a a}_{i i}^{' '} cos cos α α - - {b b}_{i i}^{' '} sin sin α α)) + + {I I}_{y the y,, i i} (({a a}_{i i}^{' '} sin sin α α \\ - - {b b}_{i i}^{' '} cos cos α α))]] \end{matrix};;$

${f f}_{3434} ((M m)) = = {Σ Σ}_{i i = = 11}^{N N} [[{I I}_{x x,, i i} ((- - {a a}_{i i}^{' '} s the s i i n no α α + + {b b}_{i i}^{' '} c c o o s the s α α)) + + {I I}_{y the y,, i i} ((- - {a a}_{i i}^{' '} c c o o s the s α α - - {b b}_{i i}^{' '} s the s i i n no α α))]];;$

f₄₄(M)＝N；f ₄₄ (M)=N;

a_i′＝a_i-a₁,b_i′＝b_i-b₁,I_t,i＝I(a_i,b_i)-I_i；a _i '=a _i -a ₁ , b _i '=b _i -b ₁ ,I _t,i =I(a _i ,b _i )-I _i ;

步骤五、计算经纬度误差的增量、航向误差的增量和地磁日变误差的增量δM：Step five, calculate the increment of latitude and longitude error, the increment of heading error and the increment δM of geomagnetic diurnal variation error:

δM＝-H×G (4)δM=-H×G (4)

步骤六、更新经纬度误差、航向误差和地磁日变误差M：Step 6. Update latitude and longitude error, heading error and geomagnetic diurnal variation error M:

M＝M+δM (5)M=M+δM (5)

步骤七、判断是否满足终止迭代条件，若满足则停止迭代并跳到步骤十，否则跳到步骤八；Step 7. Judging whether the termination iteration condition is satisfied, if so, stop the iteration and skip to step 10, otherwise skip to step 8;

终止迭代条件为①、②中任意一个或两个：①迭代次数达到预设次数；②经纬度误差的增量、航向误差的增量和地磁日变误差的增量δM的2范数小于设定值，即Termination of the iteration condition is any one or both of ① and ②: ① the number of iterations reaches the preset number; ② the increment of latitude and longitude error, the increment of heading error and the 2-norm of δM of the increment of geomagnetic diurnal variation error are less than the set value, ie

||δM||₂＜ε (6)||δM|| ₂ <ε (6)

其中ε为预先设定的迭代最小误差；Where ε is the preset iteration minimum error;

步骤八、根据更新后的M计算参数K和δK：Step 8. Calculate parameters K and δK according to the updated M:

K＝g(M) (7)K=g(M) (7)

δK＝K-G (8)δK=K-G (8)

步骤九、根据公式(9)和公式(10)更新迭代变量H和G，然后跳到步骤五，Step nine, update iteration variables H and G according to formula (9) and formula (10), then skip to step five,

H＝H+[(δM)-H(δK)](δM)^TH/[(δM)^TH(δK)] (9)H＝H+[(δM)-H(δK)](δM) ^T H/[(δM) ^T H(δK)] (9)

G＝K (10)G=K (10)

步骤十、根据迭代计算得到经纬度误差、航向误差和地磁日变误差Δx、Δy、α和δ，将所得结果代入匹配轨迹与参考轨迹的关系方程(11)即得匹配轨迹；Step ten, obtain latitude and longitude error, heading error and geomagnetic diurnal variation error Δx, Δy, α and δ according to iterative calculation, substitute the obtained result into the relationship equation (11) between the matching trajectory and the reference trajectory to obtain the matching trajectory;

$[\begin{matrix} {u u}_{i i} \\ {v v}_{i i} \end{matrix}] = = [\begin{matrix} c c o o s the s α α & s the s i i n no α α \\ - - s the s i i n no α α & c c o o s the s α α \end{matrix}] [\begin{matrix} {a a}_{i i} - - {a a}_{11} \\ {b b}_{i i} - - {b b}_{11} \end{matrix}] + + [\begin{matrix} {a a}_{11} \\ {b b}_{11} \end{matrix}] + + [\begin{matrix} Δ Δ x x \\ Δ Δ y the y \end{matrix}] - - - - - - ((1111))$

其中u_i为第i时刻匹配结果的位置经度，v_i为第i时刻匹配结果的位置纬度。Where u _i is the location longitude of the matching result at the i-th moment, and v _i is the location latitude of the matching result at the i-th moment.

本发明以惯导系统输出的参考轨迹、计算机中的地磁图、地磁传感器实测值作为输入，用地磁日变影响下的惯性地磁匹配定位方法解算出参考轨迹初始经纬度误差、航向误差，将解算得到的经纬度误差、航向误差代入匹配轨迹与参考轨迹的关系方程即得到匹配结果。本发明所提出的地磁日变影响下的惯性地磁匹配定位方法具有较高的定位精度，而且重复迭代后，定位精度得到了进一步提高。The present invention uses the reference track output by the inertial navigation system, the geomagnetic map in the computer, and the measured value of the geomagnetic sensor as input, and uses the inertial geomagnetic matching positioning method under the influence of the daily change of the geomagnetic field to solve the initial latitude and longitude error and heading error of the reference track, and solve the calculation The obtained latitude and longitude error and heading error are substituted into the relationship equation between the matching trajectory and the reference trajectory to obtain the matching result. The inertial geomagnetic matching positioning method under the influence of the diurnal variation of the geomagnetic field proposed by the present invention has high positioning accuracy, and after repeated iterations, the positioning accuracy is further improved.

附图说明Description of drawings

图1为以某实验跑车为例实施本发明方法的流程图。Fig. 1 is a flow chart of implementing the method of the present invention by taking a certain experimental sports car as an example.

图2为经度误差曲线图。Figure 2 is a graph of the longitude error.

图3为纬度误差曲线图。Figure 3 is a graph of latitude error.

具体实施方式detailed description

下面将对本发明做进一步的详细说明：本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式，但本发明的保护范围不限于下述实施例。The present invention will be described in further detail below: the present embodiment is implemented on the premise of the technical solution of the present invention, and detailed implementation is provided, but the protection scope of the present invention is not limited to the following examples.

本实施例所涉及的一种地磁日变影响下的惯性地磁匹配定位方法，包括以下步骤：An inertial geomagnetic matching and positioning method under the influence of geomagnetic diurnal variation involved in this embodiment includes the following steps:

M＝[Δx Δy α δ]^T＝[0 0 0 0]^T M＝[Δx Δy α δ] ^T ＝[0 0 0 0] ^T

G＝g(M) (1)G=g(M) (1)

F＝f(M) (2)F＝f(M) (2)

H＝F^-1 (3)H＝F ^-1 (3)

其中：in:

${f f}_{1111} ((M m)) = = {Σ Σ}_{i i = = 11}^{N N} {I I}_{x x,, i i}^{22};;$

${f f}_{1414} ((M m)) = = {Σ Σ}_{i i = = 11}^{N N} {I I}_{x x,, i i};;$

f₄₄(M)＝N；f ₄₄ (M)=N;

δM＝-H×G (4)δM=-H×G (4)

M＝M+δM (5)M=M+δM (5)

||δM||₂＜ε (6)||δM|| ₂ <ε (6)

K＝g(M) (7)K=g(M) (7)

δK＝K-G (8)δK=K-G (8)

G＝K (10)G=K (10)

前文所述地磁场强度可以选择地磁场总强度、地磁异常场总强度或者地磁场总强度在地理坐标系下某一方向的分量；当地磁场强度为地磁场总强度时，磁强计采用标量磁强计或者三轴矢量磁强计，从所述磁强计直接获得地磁场总强度，作为地磁场强度的测量值I_i；相应地，预先存储的地磁场强度和地磁场强度的梯度应为地磁场总强度和地磁场总强度的梯度。当地磁场强度为地磁异常场总强度时，磁强计采用标量磁强计或者三轴矢量磁强计，从磁强计直接获得地磁场总强度，并根据地球磁场模型计算出地磁异常场总强度，作为地磁场强度的测量值I_i；相应地，预先存储的地磁场强度和地磁场强度的梯度应为地磁异常场总强度和地磁异常场总强度的梯度。当地磁场强度为地磁场总强度在地理坐标系下某一方向的分量时，磁强计采用三轴矢量磁强计，依据磁强计的三轴矢量测量值以及载体的姿态，计算出地磁场总强度在地理坐标系下该方向的分量，作为地磁场强度的测量值I_i；相应地，预先存储的地磁场强度和地磁场强度的梯度应为地磁场总强度在地理坐标系下该方向的分量及该分量的梯度。The strength of the geomagnetic field mentioned above can be selected from the total strength of the geomagnetic field, the total strength of the geomagnetic anomaly field, or the component of the total strength of the geomagnetic field in a certain direction under the geographic coordinate system; when the local magnetic field strength is the total strength of the geomagnetic field, the magnetometer uses scalar magnetic A strong meter or a three-axis vector magnetometer, directly obtain the total strength of the geomagnetic field from the magnetometer, as the measured value I _i of the geomagnetic field strength; correspondingly, the gradient of the pre-stored geomagnetic field strength and the geomagnetic field strength should be The total strength of the Earth's magnetic field and the gradient of the total Earth's magnetic field strength. When the local magnetic field strength is the total strength of the geomagnetic anomaly field, the magnetometer adopts a scalar magnetometer or a three-axis vector magnetometer to directly obtain the total geomagnetic field strength from the magnetometer, and calculate the total geomagnetic anomaly field strength according to the earth magnetic field model , as the measured value I _i of the geomagnetic field strength; correspondingly, the pre-stored geomagnetic field strength and the gradient of the geomagnetic field strength should be the total strength of the geomagnetic anomaly field and the gradient of the total geomagnetic anomaly field strength. When the local magnetic field strength is the component of the total strength of the geomagnetic field in a certain direction under the geographic coordinate system, the magnetometer uses a three-axis vector magnetometer, and calculates the geomagnetic field according to the three-axis vector measurement value of the magnetometer and the attitude of the carrier. The component of the total strength in this direction under the geographic coordinate system is taken as the measured value I _i of the geomagnetic field strength; correspondingly, the pre-stored geomagnetic field strength and the gradient of the geomagnetic field strength should be the total strength of the geomagnetic field in this direction under the geographic coordinate system The component of and the gradient of the component.

地磁日变影响下的惯性地磁匹配定位方法还可以包括二次迭代步骤，所述二次迭代步骤的内容为：第二次执行所述步骤一到十，与第一次执行(即一次迭代)有如下两点不同：①用第一次执行所得到的匹配定位结果u_i和v_i替代所述步骤一中从惯导系统读取的a_i和b_i(如式(12)所示)，即将一次迭代得到的修正后的待匹配点位置值作为二次迭代的初始值；②用第一次执行所得到的地磁日变误差δ作为所述步骤三中的初始地磁日变误差(如式(13)所示)。由于二次迭代用的初值是经过第一次迭代获得的已经消除了大部分初始位置误差、大部分初始航向误差的值，再次通过迭代计算，其结果比第一次迭代更接近真实轨迹，使匹配定位结果精度得到进一步提高。The inertial geomagnetic matching positioning method under the influence of geomagnetic diurnal variation can also include a second iteration step, the content of the second iteration step is: the second execution of the steps one to ten, and the first execution (i.e. an iteration) There are two differences as follows: ① Replace the a _i and b _i read from the inertial navigation system in the first step with the matching positioning results u _i and v _i obtained in the first execution (as shown in formula (12)) , the corrected position value of the point to be matched obtained in the first iteration is taken as the initial value of the second iteration; ② the geomagnetic diurnal variation error δ obtained in the first execution is used as the initial geomagnetic diurnal variation error in the step 3 (such as shown in formula (13)). Since the initial value used for the second iteration is the value obtained through the first iteration that has eliminated most of the initial position error and most of the initial heading error, and through iterative calculation again, the result is closer to the real trajectory than the first iteration. The accuracy of matching positioning results is further improved.

$\{\begin{matrix} {a a}_{i i} = = {u u}_{i i} \\ {b b}_{i i} = = {v v}_{i i} \end{matrix} - - - - - - ((1212))$

M＝[0 0 0 δ]^T (13)M=[0 0 0 δ] ^T (13)

依照图1，以某实验跑车为例实施本发明地磁日变影响下的惯性地磁匹配定位方法的过程如下：According to Fig. 1, taking an experimental sports car as an example, the process of implementing the inertial geomagnetic matching positioning method under the influence of geomagnetic diurnal variation of the present invention is as follows:

实验条件：地磁场选用地磁异常场总强度。选用质子磁力仪实时测量磁场信息，质子磁力仪的主要性能指标如下：分辨率：0.01nT，精度：±0.2nT。惯导系统的主要性能指标如下：陀螺零偏不稳定性：0.01°/h，陀螺随机游走：加速度计零偏不稳定性：80μg，加速度计随机游走： Experimental conditions: The total strength of the geomagnetic anomaly field is selected for the geomagnetic field. The proton magnetometer is selected to measure the magnetic field information in real time. The main performance indicators of the proton magnetometer are as follows: resolution: 0.01nT, accuracy: ±0.2nT. The main performance indicators of the inertial navigation system are as follows: gyro bias instability: 0.01°/h, gyro random walk: Accelerometer bias instability: 80μg, accelerometer random walk:

将经度范围126°到129°和纬度范围44°到46°区间的地磁异常场总强度数据存入机载计算机，采用前向差分方法计算地磁异常场总强度的梯度信息并存入机载计算机；选取待匹配点个数为7个，即N＝7。Store the total strength data of the geomagnetic anomaly field in the longitude range from 126° to 129° and the latitude range from 44° to 46° into the onboard computer, and use the forward difference method to calculate the gradient information of the total strength of the geomagnetic anomaly field and store it in the onboard computer ;Select the number of points to be matched as 7, that is, N=7.

采用地磁日变影响下的惯性地磁匹配定位方法步骤如下：The steps of the inertial geomagnetic matching positioning method under the influence of geomagnetic diurnal variation are as follows:

步骤一：从惯性导航系统读取当前时刻以及前6个时刻的待匹配点的位置测量值a_i和b_i，如表1所示；根据质子磁力仪的测量值和地球磁场模型，得到当前时刻以及前6个时刻的磁场强度测量信息I_i如表2所示：Step 1: Read the position measurements a _i and b _i of the points to be matched at the current moment and the previous 6 moments from the inertial navigation system, as shown in Table 1; according to the measured values of the proton magnetometer and the earth's magnetic field model, the current The time and the magnetic field strength measurement information I _i of the first 6 moments are shown in Table 2:

表1惯性导航系统测量位置坐标Table 1 Inertial navigation system measurement position coordinates

ii 经度a_i(°)Longitude a _i (°) 纬度b_i(°)Latitude b _i (°) 11 127.1558127.1558 44.398044.3980 22 127.1988127.1988 44.415444.4154 33 127.2619127.2619 44.445644.4456 44 127.3237127.3237 44.474644.4746 55 127.3836127.3836 44.503344.5033 66 127.4410127.4410 44.532444.5324 77 127.4946127.4946 44.563044.5630

表2磁场强度测量信息Table 2 Magnetic Field Strength Measurement Information

ii 磁场强度I_i(nT)Magnetic field strength I _i (nT) 11 134.23134.23 22 142.36142.36 33 157.68157.68 44 148.19148.19 55 142.62142.62 66 146.95146.95 77 162.27162.27

步骤二：根据惯性导航系统指示的7个位置，分别从预先存储的地磁数据库中读取该位置的地磁场强度信息I(a_i,b_i)，以及该位置的地磁场强度的梯度信息I_x,i和I_y,i如表3所示：Step 2: According to the 7 positions indicated by the inertial navigation system, read the geomagnetic field strength information I(a _i , _bi ) of the position and the gradient information I of the geomagnetic field strength of the position from the pre-stored geomagnetic database respectively _{x, i} and I _{y, i} are shown in Table 3:

表3地磁数据库中的地磁场强度和梯度信息Table 3 Geomagnetic field strength and gradient information in the geomagnetic database

步骤三：初始化经纬度误差、航向误差以及地磁日变误差：Step 3: Initialize latitude and longitude error, heading error and geomagnetic diurnal variation error:

M＝[0 0 0 0]^T M=[0 0 0 0] ^T

步骤四：根据公式(1)、(2)和(3)计算迭代参数G、F和HStep 4: Calculate iteration parameters G, F and H according to formulas (1), (2) and (3)

$G G = = [\begin{matrix} - - 0.7557 0.7557 \\ - - 4.6567 4.6567 \\ 0.0495 0.0495 \\ - - 0.0116 0.0116 \end{matrix}] \times \times 1010^{44}$

$F f = = [\begin{matrix} 0.1878 0.1878 & 0.0494 0.0494 & 0.0237 0.0237 & 0.0008 0.0008 \\ 0.0494 0.0494 & 1.0446 1.0446 & - - 0.0328 0.0328 & 0.0014 0.0014 \\ 0.0237 0.0237 & - - 0.0328 0.0328 & 0.0092 0.0092 & 0.0001 0.0001 \\ 0.0008 0.0008 & 0.0014 0.0014 & 0.0001 0.0001 & 0.0001 0.0001 \end{matrix}] \times \times 1010^{66}$

$H h = = [\begin{matrix} 0.0161 0.0161 & 0.0013 0.0013 & 0.0201 0.0201 & 1.9347 1.9347 \\ 0.0013 0.0013 & 0.0023 0.0023 & 0.0112 0.0112 & 0.7456 0.7456 \\ 0.0201 0.0201 & 0.0112 0.0112 & 0.2207 0.2207 & 2.4284 2.4284 \\ 1.9347 1.9347 & 0.7456 0.7456 & 2.4284 2.4284 & 552.5414 552.5414 \end{matrix}] \times \times 1010^{- - 33}$

步骤五至步骤九：选取ε＝10^-6，预设迭代次数为50次。执行公式(4)、公式(5)、公式(7)至公式(10)，实施迭代算法。并根据公式(6)判断迭代终止条件，可知当迭代次数为7时迭代终止。迭代计算得到的经纬度误差、航向误差以及地磁日变误差M为：Step 5 to Step 9: select ε=10 ^-6 , and the default number of iterations is 50. Execute formula (4), formula (5), formula (7) to formula (10) to implement the iterative algorithm. And judge the iteration termination condition according to the formula (6), it can be known that the iteration terminates when the number of iterations is 7. The longitude and latitude error, heading error and geomagnetic diurnal variation error M obtained by iterative calculation are:

M＝[-0.0325 0.0225 -0.0382 16.2869]^T M＝[-0.0325 0.0225 -0.0382 16.2869] ^T

步骤十：将算得的误差值代入匹配轨迹与参考轨迹关系方程(11)，算得匹配后的经度、纬度值如表4所示：Step 10: Substituting the calculated error value into the relationship equation (11) between the matching track and the reference track, and calculating the matched longitude and latitude values as shown in Table 4:

表4匹配结果位置坐标Table 4 Location coordinates of matching results

ii 经度u_i(°)Longitude u _i (°) 纬度v_i(°)Latitude v _i (°) 11 127.1233127.1233 44.420544.4205 22 127.1656127.1656 44.439544.4395 33 127.2275127.2275 44.472144.4721 44 127.2882127.2882 44.503544.5035 55 127.3469127.3469 44.534444.5344 66 127.4032127.4032 44.565744.5657 77 127.4556127.4556 44.598344.5983

上述迭代计算在计算机中耗时为15.5毫秒。The above iterative calculation takes 15.5 milliseconds in the computer.

采用二次迭代来提高匹配定位方法的精度。二次迭代时，将步骤一中从惯导系统读取的待匹配点的位置测量值a_i和b_i由第一次迭代所得到的匹配定位结果u_i和v_i来替代，并将步骤三中初始地磁日变误差由第一次迭代所得到的地磁日变误差δ来替代。二次迭代步骤如下：A second iteration is used to improve the accuracy of the matching localization method. In the second iteration, the position measurement values a _i and b _i of the points to be matched read from the inertial navigation system in step 1 are replaced by the matching positioning results u _i and v _i obtained in the first iteration, and the step The initial geomagnetic diurnal variation error in the third step is replaced by the geomagnetic diurnal variation error δ obtained from the first iteration. The second iteration steps are as follows:

步骤一：用第一次的迭代结果(见表4)替代从惯导系统读取的7个时刻的待匹配点的位置测量值a_i和b_i，如表5所示；根据质子磁力仪的测量值和地球磁场模型，得到当前时刻以及前6个时刻的磁场强度测量信息I_i，如表2所示：Step 1: Use the first iteration result (see Table 4) to replace the position measurement values a _i and b _i of the points to be matched at seven moments read from the inertial navigation system, as shown in Table 5; according to the proton magnetometer The measured value and the earth's magnetic field model, the magnetic field strength measurement information I _i of the current moment and the previous 6 moments are obtained, as shown in Table 2:

表5二次迭代时的待匹配点位置坐标初值Table 5 The initial value of the position coordinates of the points to be matched during the second iteration

ii 经度a_i(°)Longitude a _i (°) 纬度b_i(°)Latitude b _i (°) 11 127.1233127.1233 44.420544.4205 22 127.1656127.1656 44.439544.4395 33 127.2275127.2275 44.472144.4721 44 127.2882127.2882 44.503544.5035 55 127.3469127.3469 44.534444.5344 66 127.4032127.4032 44.565744.5657 77 127.4556127.4556 44.598344.5983

步骤二：根据表5中的7个位置，从预先存储的地磁数据库中分别读取这7个位置的地磁场强度I(a_i,b_i)，以及该位置的地磁场强度的梯度I_x,i和I_y,i，如表6所示：Step 2: According to the 7 positions in Table 5, respectively read the geomagnetic field strength I(a _i , _bi ) of these 7 positions from the pre-stored geomagnetic database, and the gradient I _x of the geomagnetic field strength at this position _,i and I _y,i , as shown in Table 6:

表6地磁数据库中的地磁场强度和梯度信息Table 6 Geomagnetic field strength and gradient information in geomagnetic database

步骤三：用第一次迭代所得到的地磁日变误差δ来初始化本次地磁日变误差，并初始化经纬度误差以及航向误差：Step 3: Use the geomagnetic daily variation error δ obtained in the first iteration to initialize the current geomagnetic daily variation error, and initialize the latitude and longitude error and heading error:

M＝[0 0 0 16.2869]^T M＝[0 0 0 16.2869] ^T

$G G = = [\begin{matrix} - - 4.0307 4.0307 \\ 0.4355 0.4355 \\ - - 1.1574 1.1574 \\ - - 0.0353 0.0353 \end{matrix}] \times \times 1010^{33}$

$F f = = [\begin{matrix} 2.8352 2.8352 & - - 0.0371 0.0371 & 0.5985 0.5985 & 0.0102 0.0102 \\ - - 0.0371 0.0371 & 5.3856 5.3856 & - - 0.6101 0.6101 & - - 0.0016 0.0016 \\ 0.5985 0.5985 & - - 0.6101 0.6101 & 0.2257 0.2257 & 0.003 0.003 \\ 0.0102 0.0102 & - - 0.0016 0.0016 & 0.003 0.003 & 0.0001 0.0001 \end{matrix}] \times \times 1010^{55}$

$H h = = [\begin{matrix} 0.0174 0.0174 & - - 0.0078 0.0078 & - - 0.0709 0.0709 & 0.2682 0.2682 \\ - - 0.0078 0.0078 & 0.0076 0.0076 & 0.0540 0.0540 & - - 0.9626 0.9626 \\ - - 0.0709 0.0709 & 0.0540 0.0540 & 0.5067 0.5067 & - - 9.7977 9.7977 \\ 0.2682 0.2682 & - - 0.9626 0.9626 & - - 9.7977 9.7977 & 495.3714 495.3714 \end{matrix}] \times \times 1010^{- - 33}$

M＝[0.0055 -0.0086 -0.0903 24.2326]^T M＝[0.0055 -0.0086 -0.0903 24.2326] ^T

步骤十：将算得的误差值代入匹配轨迹与参考轨迹关系方程(11)，算得匹配后的经度、纬度值如表7所示：Step 10: Substituting the calculated error value into the relational equation (11) between the matching trajectory and the reference trajectory, the calculated longitude and latitude values after matching are shown in Table 7:

表7二次迭代匹配结果位置坐标Table 7 Position coordinates of the second iteration matching result

ii 经度u_i(°)Longitude u _i (°) 纬度v_i(°)Latitude v _i (°) 11 127.1288127.1288 44.411944.4119 22 127.1693127.1693 44.434744.4347 33 127.2280127.2280 44.472744.4727 44 127.2855127.2855 44.509444.5094 55 127.3413127.3413 44.545544.5455 66 127.3945127.3945 44.581844.5818 77 127.4437127.4437 44.619044.6190

两次迭代计算在计算机中总耗时为17.1毫秒。The total time spent in the computer for two iterations is 17.1 milliseconds.

为了说明校正地磁日变后的效果，我们还得到不校正地磁日变时的匹配结果位置坐标，如表8所示：In order to illustrate the effect of correcting the geomagnetic diurnal variation, we also obtained the position coordinates of the matching results without correcting the geomagnetic diurnal variation, as shown in Table 8:

表8无地磁日变校正的匹配位置坐标Table 8 Matching position coordinates without geomagnetic diurnal correction

ii 经度u_i′(°)Longitude u _i ′(°) 纬度v_i′(°)Latitude v _i ′(°) 11 127.1841127.1841 44.437344.4373 22 127.2299127.2299 44.444544.4445 33 127.2982127.2982 44.459644.4596 44 127.3650127.3650 44.473944.4739 55 127.4298127.4298 44.488344.4883 66 127.4923127.4923 44.503644.5036 77 127.5515127.5515 44.521344.5213

为了验证实验结果，在实验跑车上安装GPS导航定位系统，从而可以得到7个时刻的真实位置坐标，如表9所示：In order to verify the experimental results, a GPS navigation and positioning system is installed on the experimental sports car, so that the real position coordinates at seven moments can be obtained, as shown in Table 9:

表9跑车真实位置坐标Table 9 The real position coordinates of the sports car

ii 经度(°)Longitude (°) 纬度(°)Latitude (°) 11 127.1372127.1372 44.409444.4094 22 127.1680127.1680 44.431444.4314 33 127.2232127.2232 44.471044.4710 44 127.2786127.2786 44.510644.5106 55 127.3339127.3339 44.550244.5502 66 127.3893127.3893 44.589844.5898 77 127.4447127.4447 44.629444.6294

根据GPS导航定位系统的定位结果(表9)，可以绘出7个时刻的无日变校正的定位误差曲线、有日变校正的第一次迭代结果的定位误差曲线和有日变校正的第二次迭代结果的定位误差曲线，如图2和图3所示，其中图2给出了经度误差曲线，图3给出了纬度误差曲线。根据图2和图3可以看出，所提出的地磁日变影响下的惯性地磁匹配定位方法具有较高的定位精度，而且重复迭代后，定位精度会进一步提高。According to the positioning results of the GPS navigation and positioning system (Table 9), the positioning error curve without diurnal variation correction, the positioning error curve of the first iteration result with diurnal variation correction, and the first iteration with diurnal variation correction can be drawn at seven times. The positioning error curve of the second iteration result is shown in Figure 2 and Figure 3, where Figure 2 shows the longitude error curve, and Figure 3 shows the latitude error curve. According to Fig. 2 and Fig. 3, it can be seen that the proposed inertial geomagnetic matching positioning method under the influence of geomagnetic diurnal variation has high positioning accuracy, and after repeated iterations, the positioning accuracy will be further improved.

以上所述，仅为本发明较佳的具体实施方式，这些具体实施方式都是基于本发明整体构思下的不同实现方式，而且本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求书的保护范围为准。The above are only preferred specific implementations of the present invention. These specific implementations are all based on different implementations under the overall concept of the present invention, and the scope of protection of the present invention is not limited thereto. Anyone familiar with the technical field Within the technical scope disclosed in the present invention, any changes or substitutions that can be easily conceived by a skilled person shall fall within the protection scope of the present invention. Therefore, the protection scope of the present invention should be determined by the protection scope of the claims.

Claims

1. the inertia earth magnetism matching locating method under the influence of a geomagnetic diurnal change, it is characterised in that

Step one, read current time and the position measurements a of the point to be matched in front N-1 moment from inertial navigation system_iAnd b_i, Wherein a_iRepresent longitude, b_iRepresenting latitude, subscript i represents the most in the same time, i=1 ... N, N are integer and N ＞ 2, and i is that 1 expression is current In the moment, obtain current time and the measured value I of the geomagnetic field intensity in front N-1 moment by geomagnetic sensor_i；

Step 2, the position of the point N number of to be matched indicated according to inertial navigation system, read from the geomagnetic database prestored respectively Take the reference value I (a of corresponding geomagnetic field intensity_i,b_i), reference value I of the gradient of geomagnetic field intensity_x,iAnd I_y,i, wherein I_x,iTable Show the geomagnetic field intensity gradient in longitudinal value on i-th position, I_y,iRepresent that geomagnetic field intensity is latitudinal Gradient value on i-th position；

Step 3, introduce and initialize longitude and latitude error, course error and geomagnetic diurnal change error:

M=[Δ x Δ y α δ]^T=[0 00 0]^T

Step 4, calculate iterative parameter G, F and H according to formula (1), formula (2) and formula (3):

G=g (M) (1)

F=f (M) (2)

H=F^-1 (3)

Wherein:

G (M)=[g₁(M) g₂(M) g₃(M) g₄(M)]^T；

f (M) = [\begin{matrix} f_{11} (M) & f_{12} (M) & f_{13} (M) & f_{14} (M) \\ f_{12} (M) & f_{22} (M) & f_{23} (M) & f_{24} (M) \\ f_{13} (M) & f_{23} (M) & f_{33} (M) & f_{34} (M) \\ f_{14} (M) & f_{24} (M) & f_{34} (M) & f_{44} (M) \end{matrix}];

\begin{matrix} g_{1} (M) = Σ_{i = 1}^{N} I_{x, i} [I_{x, i} (a_{i}^{'} \cos α + b_{i}^{'} \sin α - a_{i}^{'} + Δ x) \\ + I_{y, i} (- a_{i}^{'} \sin α + b_{i}^{'} \cos α - b_{i}^{'} + Δ y) + I_{t, i} + δ] \end{matrix};

\begin{matrix} g_{2} (M) = Σ_{i = 1}^{N} I_{y, i} [I_{x, i} (a_{i}^{'} \cos α + b_{i}^{'} \sin α - a_{i}^{'} + Δ x) \\ + I_{y, i} (- a_{i}^{'} \sin α + b_{i}^{'} \cos α - b_{i}^{'} + Δ y) + I_{t, i} + δ] \end{matrix};

\begin{matrix} g_{3} (M) = Σ_{i = 1}^{N} [I_{x, i} (a_{i}^{'} \cos α + b_{i}^{'} \sin α - a_{i}^{'} + Δ x) + I_{y, i} (- a_{i}^{'} \sin α \\ + b_{i}^{'} \cos α - b_{i}^{'} + Δ y) + I_{t, i} + δ] \times [I_{x, i} (- a_{i}^{'} \sin α + b_{i}^{'} \cos α) \\ + I_{y, i} (- a_{i}^{'} \cos α - b_{i}^{'} \sin α)] \end{matrix};

\begin{matrix} g_{4} (M) = Σ_{i = 1}^{N} [I_{x, i} (a_{i}^{'} \cos α + b_{i}^{'} \sin α - a_{i}^{'} + Δ x) \\ + I_{y, i} (- a_{i}^{'} \sin α + b_{i}^{'} \cos α - b_{i}^{'} + Δ y) + I_{t, i} + δ] \end{matrix};

f_{11} (M) = Σ_{i = 1}^{N} I_{x, i}^{2};

f_{12} (M) = Σ_{i = 1}^{N} I_{x, i} I_{y, i};

f_{13} (M) = Σ_{i = 1}^{N} I_{x, i} [I_{x, i} (- a_{i}^{'} \sin α + b_{i}^{'} \cos α) + I_{y, i} (- a_{i}^{'} \cos α - b_{i}^{'} \sin α)];

f_{14} (M) = Σ_{i = 1}^{N} I_{x, i};

f_{22} (M) = Σ_{i = 1}^{N} I_{y, i}^{2};

f_{23} (M) = Σ_{i = 1}^{N} I_{y, i} [I_{x, i} (- a_{i}^{'} \sin α + b_{i}^{'} \cos α) + I_{y, i} (- a_{i}^{'} \cos α - b_{i}^{'} \sin α)];

f_{24} (M) = Σ_{i = 1}^{N} I_{y, i};

\begin{matrix} f_{33} (M) = Σ_{i = 1}^{N} {[I_{x, i} (- a_{i}^{'} \sin α + b_{i}^{'} \cos α) + I_{y, i} (- a_{i}^{'} \cos α - b_{i}^{'} \sin α)]}^{2} \\ + [I_{x, i} (a_{i}^{'} \cos α + b_{i}^{'} \sin α - a_{i}^{'} + Δ x) + I_{y, i} (- a_{i}^{'} \sin α + b_{i}^{'} \cos α \\ - b_{i}^{'} + Δ y) + I_{t, i} + δ] \times [I_{x, i} (- a_{i}^{'} \cos α - b_{i}^{'} \sin α) + I_{y, i} (a_{i}^{'} \sin α \\ - b_{i}^{'} \cos α)] \end{matrix};

f_{34} (M) = Σ_{i = 1}^{N} [I_{x, i} (- a_{i}^{'} s i n α + b_{i}^{'} c o s α) + I_{y, i} (- a_{i}^{'} c o s α - b_{i}^{'} s i n α)];

f₄₄(M)=N；

a′_i=a_i-a₁, b '_i=b_i-b₁, I_t,i=I (a_i,b_i)-I_i；

Step 5, the calculating increment of longitude and latitude error, the increment of course error and the increment δ M of geomagnetic diurnal change error:

δ M=-H × G (4)

Step 6, renewal longitude and latitude error, course error and geomagnetic diurnal change error M:

M=M+ δ M (5)

Step 7, judge whether to meet and terminate iterated conditional, if meeting, stopping iteration and also jumping to step 10, otherwise jumping to step Eight；

Terminate iterated conditional for 1., 2. in any one or two: 1. iterations reaches preset times；2. longitude and latitude error 2 norms of the increment δ M of increment, the increment of course error and geomagnetic diurnal change error are less than setting value, i.e.

||δM||₂＜ ε (6)

Wherein ε is iteration minimum error set in advance；

Step 8, calculate parameter K and δ K according to the M after updating:

K=g (M) (7)

δ K=K-G (8)

Step 9, update iteration variable H and G according to formula (9) and formula (10), then jump to step 5,

H=H+ [(δ M)-H (δ K)] (δ M)^TH/[(δM)^TH(δK)] (9)

G=K (10)

Step 10, the longitude and latitude error that obtains according to iterative computation, course error and geomagnetic diurnal change error delta x, Δ y, α and δ, by institute Result substitutes into coupling track and i.e. obtains with the relation equation (11) of reference locus and mate track；

[\begin{matrix} u_{i} \\ v_{i} \end{matrix}] = [\begin{matrix} c o s α & s i n α \\ - s i n α & c o s α \end{matrix}] [\begin{matrix} a_{i} - a_{1} \\ b_{i} - b_{1} \end{matrix}] + [\begin{matrix} a_{1} \\ b_{1} \end{matrix}] + [\begin{matrix} Δ x \\ Δ y \end{matrix}] - - - (11)

Wherein u_iIt is the position longitude of the i-th moment matching result, v_iIt it is the position latitude of the i-th moment matching result.