CN104061945B - Deviation of plumb line dynamic measurement device based on INS and GPS combination and method - Google Patents

Abstract

The invention discloses a kind of deviation of plumb line dynamic measurement device based on INS and GPS combination and method.First constructing INS/GPS attitude measurement subsystem and LGU/GPS attitude measurement subsystem, when starting to measure, LGU/GPS attitude measurement subsystem is entered to initialize by the attitude matrix utilizing INS/GPS to export；Hereafter, utilize LGU/GPS and INS/GPS to carry out attitude algorithm respectively, and calculate the difference of the attitude angle of LGU/GPS Yu INS/GPS output；By setting up observational equation and the state equation that the deviation of plumb line is measured, with the difference of attitude angle of INS/GPS Yu LGU/GPS output as observed quantity, under the auxiliary of global Gravity Models, the method for Kalman filter is utilized to obtain the optimal estimation value of the deviation of plumb line.The present invention can effectively reduce complexity and the cost of system, improves precision and the reliability of measurement, does not relies on measurement environment, workable and convenient to carry out.

Description

Vertical line deviation dynamic measurement device and method based on INS and GPS combination

Technical Field

The invention belongs to the field of surveying and mapping, and particularly relates to a vertical line deviation dynamic measurement method based on the combination of an Inertial Navigation System (INS) and a Global Positioning System (GPS).

Background

The vector gravity measurement method based on the combination of the inertial navigation system and the global positioning system is an effective means for acquiring gravity field information, has the characteristics of high precision, high resolution, high measurement efficiency and the like, and is widely applied to marine and aviation gravity measurement. The measurement of the vertical component of the gravity vector is mature at present and widely used in commercial measurement, however, the measurement of the horizontal component of the gravity vector (i.e. the deviation of the vertical line) is still a technical problem which is difficult to solve.

The warrior and German force of national defense science and technology university, in the invention patent application with the application number of "CN 201310730904.5" (publication number: CN103674030A), provides a device and a method for measuring the vertical line deviation based on astronomical attitude reference maintenance, the adopted technical scheme is shown in figure 1, and a laser gyro assembly (LGU) in an INS/GPS attitude measurement system is combined with a GPS to construct an LGU/GPS attitude measurement subsystem; firstly, starting an INS/GPS attitude measurement system, performing initial alignment and outputting a system attitude; then starting the star sensor and utilizing the attitude output of the star sensor to initialize the LGU/GPS attitude measurement subsystem, and then enabling the LGU/GPS attitude measurement subsystem to autonomously realize attitude measurement; calculating the difference between the posture output by the LGU/GPS posture measurement subsystem and the posture output by the INS/GPS posture measurement system, and further calculating the deviation of the vertical line; and finally, correcting jump errors in the vertical deviation measurement values, and correcting low-frequency errors in the vertical deviation measurement values by using global gravity model data. The method can realize dynamic measurement of the vertical line deviation, has strong robustness, can effectively inhibit errors of inertial devices, does not depend on a differential GPS, and has the advantages of wide application range and the like.

However, the above measurement method has certain limitations, which are mainly reflected in the following three aspects: firstly, the star sensor is required to provide a high-precision astronomical attitude reference in the measuring method, so that the cost of measuring equipment and the complexity of a structure are greatly increased; secondly, the use of the star sensor depends on weather conditions and can only be used at night, so that the flexibility of measurement implementation is severely limited; secondly, the low-frequency error compensation algorithm for the vertical line deviation provided by the invention is only an engineering experience processing method, and the algorithm is not optimal, so that the compensation precision is not high, and the reliability is poor.

How to realize the measurement of the vertical deviation under the dynamic condition, inhibit the measurement error, improve the measurement precision and reliability, and reduce or even eliminate the dependence of the measurement equipment on the environment is a technical problem which is of great concern to the technicians in the field.

Disclosure of Invention

The invention provides an improvement of the invention CN201310730904.5, and aims to provide a device and a method for realizing the dynamic measurement of vertical line deviation by using the combination of INS and GPS under the condition of not depending on astronomical attitude reference, thereby improving the measurement precision and reliability, simplifying the complexity of a measurement system, not depending on the measurement environment and facilitating the implementation.

In order to solve the technical problems, the technical scheme adopted by the invention is as follows: a vertical line deviation dynamic measuring device based on INS and GPS combination comprises an Inertial Navigation System (INS), a GPS antenna, a GPS receiver and a data processing computer, wherein the INS, the GPS antenna and the GPS receiver form an INS/GPS attitude measurement subsystem (INS/GPS for short), the INS comprises three orthogonally-mounted laser gyros which are called a laser gyro assembly (LGU), the LGU, the GPS antenna and the GPS receiver form an LGU/GPS attitude measurement subsystem (LGU/GPS for short, and the LGU/GPS is described in CN201310730904.5 of the invention), and the three laser gyros are all communicated with the GPS receiver; the INS, the GPS antenna and the GPS receiver are fixedly connected and installed on a measuring carrier, and the measuring carrier can be a measuring ship, a measuring vehicle and other carrying tools; the GPS antenna is communicated with the GPS receiver, the INS and the GPS receiver are connected with a data processing computer through data lines, the measurement data of the INS and the GPS receiver are transmitted to the data processing computer through the data lines, and the calculation of the vertical line deviation is completed in the data processing computer. The INS adopts a single-axis rotary structure disclosed in the paper of "laser gyro single-axis rotary inertial navigation system" in the 2 nd of "chinese technical science", 2010, by longxingwu et al.

The invention also provides a method for dynamically measuring the deviation of the perpendicular line by using the device, which specifically comprises the following steps:

1. and constructing an INS/GPS attitude measurement subsystem, and aligning for more than 8 hours.

An INS/GPS attitude measurement subsystem is constructed by an INS, a GPS antenna and a GPS receiver, the INS/GPS attitude measurement subsystem is started to carry out alignment for more than 8 hours, the INS/GPS attitude measurement subsystem adopts the algorithm described in the step 1 implemented by the patent CN201310730904.5 in the alignment process to realize combined attitude measurement, and the INS/GPS attitude measurement subsystem continuously outputs an attitude matrix of an INS coordinate system (a system b) relative to a computed navigation coordinate system (a system n') in the whole measurement processAnd no vertical deviation data is output in the alignment process.

2. And after the INS/GPS alignment is finished, starting the LGU/GPS attitude measurement subsystem, and initializing an attitude matrix for the LGU/GPS by using the attitude matrix output after the INS/GPS alignment in the first step.

Note t_iI is more than or equal to 1 and less than or equal to N at the moment corresponding to the ith measurement sampling point, N is the total number of sampling points in the whole measurement process, t_iThe attitude matrix output by the INS/GPS at the moment isLGU/GPS output attitude matrix ofThe superscript n denotes the real navigation coordinate system (n series). The moment of starting the LGU/GPS attitude measurement subsystem is recorded as t₀I.e. i is 0, the LGU/GPS attitude matrix is initialized by: order to

3. And the LGU/GPS attitude measurement subsystem carries out attitude updating.

The LGU/GPS attitude measurement subsystem adopts special technologyThe method described in step 5.3 in CN201310730904.5 is used for attitude updating, and simultaneously, the data processing computer stores attitude matrixes output by INS/GPS at all measurement sampling momentsAnd LGU/GPS output attitude matrix

4. Using t in the whole measurement process_iAttitude matrix output by time INS/GPSAnd LGU/GPS output attitude matrixComputing attitude matricesAnd calculates the difference between the corresponding attitude angles of its outputs.

The calculation method comprises the following steps:

$C_{n^{\′}\left(i\right)}^{n\left(i\right)}=C_{b\left(i\right)}^{n\left(i\right)}\·{\left[C_{b\left(i\right)}^{n^{\′}\left(i\right)}\right]}^T---\left(1\right)$

wherein]^TDenotes the transpose of the matrix and the sign-denotes the matrix multiplication.

For convenience of description, the n' to n attitude matrix calculated at any time t is abbreviated asThe difference between the three attitude angles output by the INS/GPS and the LGU/GPS at any time t is delta phi_E，ΔΦ_N，ΔΦ_UWhere subscripts E, N, U represent the east, north, and sky components, respectively. Delta phi_E、ΔΦ_NAnd Δ Φ_UCalculated from equation (2):

$\left\{\begin{array}{c}{\mathrm{\Δ\Φ}}_E=C_{n^{\′}}^n\left(\mathrm{2,3}\right)\\ {\mathrm{\Δ\Φ}}_N=C_{n^{\′}}^n\left(\mathrm{3,1}\right)\\ {\mathrm{\Δ\Φ}}_U=C_{n^{\′}}^n\left(\mathrm{1,2}\right)\end{array}\right.---\left(2\right)$

whereinTo representThe elements of row 2, column 3,to representThe elements of row 3, column 1,to representRow 1, column 2 elements of the matrix.

5. The difference delta phi of the three attitude angles output by the INS/GPS and the LGU/GPS_E，ΔΦ_N，ΔΦ_UFor the observed quantity, the vertical deviation is extracted by a Kalman filtering method with the assistance of a global gravity model through establishing an observation equation and a state equation of the vertical deviation measurement. The specific implementation method is as follows:

the deviation of the true east perpendicular line of the position of the measurement carrier at the moment of t is recorded as η, the deviation of the true north perpendicular line is recorded as ξ, and as shown in fig. 3, the deviation values of the east perpendicular line and the north perpendicular line calculated by the global gravity model are respectively recorded asThe global gravity model is EGM2008 global gravity model, and the calculation program and the use method of the global gravity model can be selected from "http:// earth-info.nga.mil/GandG/wgs84/gravitymod/index.htmlThe deviation errors of the east and north vertical lines calculated by the global gravity model of EGM2008 are respectively recorded as η and ξ, and the relationship of equation (3) exists:

$\left\{\begin{array}{c}\mathrm{\η}=\widehat{\mathrm{\η}}+\mathrm{\δ\η}\\ \mathrm{\ξ}=\widehat{\mathrm{\ξ}}+\mathrm{\δ\ξ}\end{array}\right.---\left(3\right)$

the attitude error angle of INS/GPS output relative to n' system is recorded as v_E，v_N，v_UThe LGU/GPS output has an attitude error angle phi with respect to the n system_E、φ_N、φ_USubscript E, N, U denotes the east, north, and sky components, respectively.

5.1 an equation of state for the vertical deviation measurement is established.

Selecting a state variable of the plumb line deviation measurement system as phi_E、φ_N、φ_U、η、ξ、_UWherein_URespectively carrying out dynamic modeling on the state variables for the equivalent day zero offset of the laser gyro, phi_E、φ_N、φ_UThe following differential equation is satisfied:

$\left\{\begin{array}{c}{\stackrel{\·}{\mathrm{\φ}}}_E=-{\mathrm{\ω}}_{\mathrm{ie}}\cos L\·{\mathrm{\φ}}_U+{\mathrm{\ω}}_{\mathrm{ie}}\sin L\·{\mathrm{\φ}}_N\\ {\stackrel{\·}{\mathrm{\φ}}}_N=-{\mathrm{\ω}}_{\mathrm{ie}}\sin L\·{\mathrm{\φ}}_E\\ {\stackrel{\·}{\mathrm{\φ}}}_U={\mathrm{\ω}}_{\mathrm{ie}}\cos L\·{\mathrm{\φ}}_E-{\mathrm{\ϵ}}_U\end{array}\right.---\left(4\right)$

wherein, ω is_ieIs the rotational angular velocity of the earth, and L is the geographical latitude of the measuring point.

_UModeling is a random constant model, and then:

${\stackrel{\·}{\mathrm{\ϵ}}}_U=0---\left(5\right)$

the statistical models of the deviation errors η, ξ of the east and north vertical lines of the EGM2008 are given by the expressions (6) and (7), respectively:

$\left[\begin{array}{c}{\stackrel{\·}{x}}_E\\ \mathrm{\δ}\stackrel{\·}{\mathrm{\η}}\end{array}\right]=\left[\begin{array}{cc}0& 1\\ -{\mathrm{\ω}}_0^2& -2\mathrm{\ζ}{\mathrm{\ω}}_0\end{array}\right]\left[\begin{array}{c}{x}_E\\ \mathrm{\δ\η}\end{array}\right]+\left[\begin{array}{c}0\\ w_E\end{array}\right]---\left(6\right)$

$\left[\begin{array}{c}{\stackrel{\·}{x}}_N\\ \mathrm{\δ}\stackrel{\·}{\mathrm{\ξ}}\end{array}\right]=\left[\begin{array}{cc}0& 1\\ -{\mathrm{\ω}}_0^2& -2\mathrm{\ζ}{\mathrm{\ω}}_0\end{array}\right]\left[\begin{array}{c}{x}_N\\ \mathrm{\δ\ξ}\end{array}\right]+\left[\begin{array}{c}0\\ w_N\end{array}\right]---\left(7\right)$

wherein x_E、x_NIs an intermediate state variable; omega₀Is the natural frequency, ω, of the statistical model₀With a fixed relationship omega to the speed of movement V of the measurement carrier₀2 pi × V/10000, zeta is the damping coefficient of the statistical model, w_EAnd w_NProcess noise, w, of the statistical model for the east and north plumb deviations, respectively_E～(0,q_E)，w_N～(0,q_N) I.e. w_EAnd w_NRespectively obey a Gaussian distribution with a mean value of 0 and a variance of q_EAnd q is_N。

The state space vector of the selected vertical deviation measurement system is as follows:

x＝[φ_E,φ_N,φ_U,_U,x_E,x_N,η,ξ]^T(8)

writing equations (4) - (7) as a unified equation of state form:

$\stackrel{\·}{x}=F\·x+w---\left(9\right)$

wherein, the matrix $F=\left[\begin{array}{cccccccc}0& {\mathrm{\ω}}_{\mathrm{ie}}\sin L& -{\mathrm{\ω}}_{\mathrm{ie}}\cos L& 0& 0& 0& 0& 0\\ -{\mathrm{\ω}}_{\mathrm{ie}}\sin L& 0& 0& 0& 0& 0& 0& 0\\ {\mathrm{\ω}}_{\mathrm{ie}}\cos L& 0& 0& -1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 1\\ 0& 0& 0& 0& -{\mathrm{\ω}}_0^2& 0& -2\mathrm{\ζ}{\mathrm{\ω}}_0& 0\\ 0& 0& 0& 0& 0& -{\mathrm{\ω}}_0^2& 0& -2\mathrm{\ζ}{\mathrm{\ω}}_0\end{array}\right],$

w＝[0,0,0,0,0,0,w_E,w_N]^T. Equation (9) is the state equation of the perpendicular deviation measurement system. The process noise w — (0, Q), i.e., w follows a gaussian distribution with a mean of 0 and a variance of Q. $Q=\left[\begin{array}{cc}{O}_{6\×6}& O_{6\×2}\\ O_{2\×6}& \mathrm{diag}[q_E,q_N]\end{array}\right],$ Wherein O is_6×6、O_2×6、O_6×2A zero matrix representing 6 rows and 6 columns, 2 rows and 6 columns, 6 rows and 2 columns, respectively, diag [ q ]_E,q_N]Is represented by q_EAnd q is_NIs a diagonal matrix of diagonal elements.

Discretizing equation of state (9) yields:

x_i＝M_i/i-1·x_i-1+w_i(10)

x_idenotes x is at t_iSampled value of time, M_i/i-1Is t_i-1To t_iState transition matrix of time, M_i/i-1≈I_8×8+Δt·F_iIn which I_8×8Is an 8-dimensional identity matrix, F_iFor the matrix F at t_iThe sampled value at the time, Δ t, is the sampling interval. w is a_iObedience mean 0, variance Q_iGaussian distribution of (Q)_i＝Δt²·Q。

5.2 establishing an observation equation of the vertical deviation measurement.

The difference delta phi between the attitude angles output by INS/GPS and LGU/GPS_E，ΔΦ_NFor the observed quantity, the following observation equation is established:

$\left\{\begin{array}{c}{\mathrm{\Δ\Φ}}_E=\mathrm{\ξ}-{\mathrm{\φ}}_E+v_E\\ {\mathrm{\Δ\Φ}}_N=-\mathrm{\η}-{\mathrm{\φ}}_N+v_N\end{array}\right.---\left(11\right)$

substituting formula (3) into formula (11) to obtain:

$\left\{\begin{array}{c}{\mathrm{\Δ\Φ}}_E-\widehat{\mathrm{\ξ}}=\mathrm{\δ\ξ}-{\mathrm{\φ}}_E+v_E\\ {\mathrm{\Δ\Φ}}_N+\widehat{\mathrm{\η}}=-\mathrm{\δ\η}-{\mathrm{\φ}}_N+v_N\end{array}\right.---\left(12\right)$

selecting a new observation vector as follows:

$y={[{\mathrm{\Δ\Φ}}_E-\widehat{\mathrm{\ξ}},{\mathrm{\Δ\Φ}}_N+\widehat{\mathrm{\η}}]}^T---\left(13\right)$

rewrite equation (12) to:

y＝H·x+v (14)

wherein the matrix $H=\left[\begin{array}{cccccccc}-1& 0& 0& 0& 0& 0& 0& 1\\ 0& -1& 0& 0& 0& 0& -1& 0\end{array}\right],$ v＝[v_E,v_N]^T。

Equation (14) is an observation equation of the perpendicular deviation measurement system. INS/GPS attitude error angle v_E、v_NObserved noise, v, of the east and north observed components, respectively_E～(0,r_E)，v_N～(0,r_N) I.e. v_EAnd v_NRespectively obey a Gaussian distribution with a mean value of 0 and a variance of r_EAnd r_N. The observation noise v — (0, R), i.e. v obeys a gaussian distribution with mean 0 and variance R, the observation noise variance matrix $r=\left[\begin{array}{cc}{r}_E& 0\\ 0& r_N\end{array}\right].$

Discretizing the observation equation to obtain:

y_i＝H_i·x_i+v_i(15)

y_i、v_i、H_iare respectively asy, v, H at t_iSampled value of time, v_iObedience mean 0, variance R_iOf Gaussian distribution of R_i＝Δt²·R，H_i＝H。

5.3 estimating the state vector by using Kalman filtering algorithm according to the state equation and observation equation of the discrete form of the perpendicular deviation measuring system described by the equations (10) and (15).

The iterative algorithm of the Kalman filter algorithm is as follows:

$\left\{\begin{array}{c}{\hat{x}}_{i/i-1}=M_{i/i-1}\·{\hat{x}}_{i-1}\\ P_{i/i-1}=M_{i/i-1}\·P_{i-1}{M}_{i/i-1}^T+Q_i\\ K_i=P_{i/i-1}\·H_i^T\·{[H_i\·P_{i/i-1}\·H_i^T+R_i]}^{-1}\\ P_i=(I-K_i\·H_i)\·P_{i/i-1}\·{[I-K_i\·H_i]}^{-1}+K_i\·P_{i/i-1}\·K_i^T\\ {\hat{x}}_i={\hat{x}}_{i/i-1}+K_i\·(y_i-H_i\·{\hat{x}}_{i/i-1})\end{array}\right.---\left(16\right)$

whereinIs x_iIs determined by the estimated value of (c),is x_iOne step of predicting value of, P_iIs x_iIs estimated by the covariance matrix, P_i/i-1Is P_iOne step of predicting value of, K_iFor filter gain]^-1The symbol represents the inversion operation of the matrix.

t_iThe estimates of the errors η, ξ for the east and north plumb line deviations at the times areFrom t_iState vector estimation at timeThe following are given:

$\left\{\begin{array}{c}\mathrm{\δ}{\widehat{\mathrm{\η}}}_i={\hat{x}}_i\left(7\right)\\ \mathrm{\δ}{\widehat{\mathrm{\ξ}}}_i={\hat{x}}_i\left(8\right)\end{array}\right.---\left(17\right)$

wherein,respectively representing state vector estimation valuesThe 7 th and 8 th elements of (1).

5.4 calculate vertical deviation measurements.

t_iTime east and north plumb deviation measurement η_i、ξ_iCalculated from equations (18) and (19), respectively:

${\mathrm{\η}}_i={\widehat{\mathrm{\η}}}_i+\mathrm{\δ}{\widehat{\mathrm{\η}}}_i---\left(18\right)$

${\mathrm{\ξ}}_i={\widehat{\mathrm{\ξ}}}_i+\mathrm{\δ}{\widehat{\mathrm{\ξ}}}_i---\left(19\right)$

whereinT respectively calculated for global gravity model_iThe time measures the value of the east and north vertical deviation at the point.

The invention has the following technical effects:

1. and the attitude matrix output by the INS/GPS is utilized to carry out attitude initialization on the LGU/GPS, so that the dependence on a star sensor is eliminated, the cost is reduced, the complexity of a system is reduced, and the limitations of a measuring environment and measuring time are eliminated.

2. The state equation and the observation equation of the vertical deviation measurement system are established, the optimal estimation of the vertical deviation is realized by a Kalman filtering method under the assistance of a global gravity model, and compared with the method of the invention CN201310730904.5, the method effectively separates the low-frequency error in the measurement result and improves the precision and the reliability of the measurement method.

3. Compared with the invention CN201310730904.5, the invention only needs to initialize the posture once for the LGU/GPS, simplifies the operation and improves the measuring efficiency.

Drawings

FIG. 1 is a schematic structural diagram of a vertical line deviation dynamic measurement device disclosed in Chinese patent application CN 201310730904.5;

in the figure: INS, GPS antenna, GPS receiver, star sensor, measuring carrier, data processing computer, etc. 1, INS, 2, GPS antenna, 3, GPS receiver, 4, 5, measuring carrier and 6

FIG. 2 is a schematic view of a vertical deviation dynamic measurement apparatus according to the present invention;

in the figure: INS, GPS antenna, GPS receiver, data processing computer, and measuring carrier

FIG. 3 is a relationship of a navigational coordinate system to vertical deviation;

fig. 4 is a flow chart of the vertical deviation dynamic measurement method of the invention.

Detailed Description

The present invention will be described in further detail with reference to the accompanying drawings and examples.

As shown in fig. 2, the vertical deviation dynamic measuring apparatus used in the present invention is composed of an INS1, a GPS antenna 2, a GPS receiver 3, a data processing computer 4, and a measurement carrier 5. Wherein the INS1 and the GPS antenna 2 are respectively fixedly connected and mounted on a measurement carrier 5, and the measurement carrier may be a measurement ship, a measurement vehicle, or other vehicles. The measurement data of the INS1 and the GPS receiver 3 are transmitted to the data processing computer 4 via data lines, and the vertical line deviation is resolved in the data processing computer 4. The INS adopts a single-axis rotary structure disclosed in the paper of "laser gyro single-axis rotary inertial navigation system" in the 2 nd of "chinese technical science", 2010, by longxingwu et al.

As shown in FIG. 3, the INS coordinate system is defined as the b-system O_b-x_by_bz_bThe east-north-sky coordinate system O-xyz is defined as the real navigational coordinate system, i.e. the n system, where the O-xy plane is parallel to the earth's reference ellipsoid. The definition of the calculated navigation coordinate system is O ' -x ' y ' z ', namely n ', and the definition of the calculated navigation coordinate system can be described with reference to pages 225 to 226 of inertial navigation system technology published by university press, high concentrationp, 2012. Due to the existence of the deviation of the vertical line, the calculated navigation coordinate system n' obtained by inertial navigation calculation has deviation from the real navigation coordinate system n.

FIG. 4 is a general flowchart of the measurement method of the present invention, which comprises the following steps:

step 1: and constructing an INS/GPS combined attitude measurement system, and aligning for more than 8 hours.

By INS1, GPS antenna 2 and in the arrangement of FIG. 2The GPS receiver 3 builds an INS/GPS attitude measurement system, firstly starts an INS/GPS attitude measurement subsystem, and carries out alignment for more than 8 hours. The INS/GPS combination adopts the patent "CN 201310730904.5" to implement the algorithm described in the step 1 to realize the attitude measurement. The INS/GPS attitude measurement subsystem continuously outputs an attitude matrix of an INS coordinate system (a system b) relative to a computed navigation coordinate system (a system n') in the whole measurement processAnd no vertical deviation data is output in the alignment process.

Step 2: and after the INS/GPS alignment is finished, starting the LGU/GPS attitude measurement subsystem, and initializing an attitude matrix for the LGU/GPS by using the attitude matrix output by the INS/GPS in the first step.

Note t_iI is more than or equal to 1 and less than or equal to N at the moment corresponding to the ith measurement sampling point, N is the total number of sampling points in the whole measurement process, t_iThe attitude matrix output by the INS/GPS at the moment isLGU/GPS output attitude matrix ofThe superscript n denotes the real navigation coordinate system (n series). The moment of starting the LGU/GPS attitude measurement subsystem is recorded as t₀The LGU/GPS attitude matrix initialization method is: order to

And step 3: and the LGU/GPS attitude measurement subsystem carries out attitude updating.

The LGU/GPS attitude measurement subsystem adopts the method described in step 5.3 in the patent CN201310730904.5 to update the attitude, and simultaneously the data processing computer stores the attitude matrix output by INS/GPS at all the measurement sampling momentsAnd LGU/GPS output attitude matrix

And 4, step 4: using t in the whole measurement process_iAttitude matrix output by time INS/GPSAnd LGU/GPS output attitude matrixComputing attitude matricesAnd calculates the difference between the corresponding attitude angles of its outputs.

The calculation method comprises the following steps:

$C_{n^{\′}\left(i\right)}^{n\left(i\right)}=C_{b\left(i\right)}^{n\left(i\right)}\·{\left[C_{b\left(i\right)}^{n^{\′}\left(i\right)}\right]}^T---\left(1\right)$

wherein]^TDenotes the transpose of the matrix and the sign-denotes the matrix multiplication.

Simplifying the attitude matrix from n' system to n system obtained by calculating at any time t asThe difference between the three attitude angles output by the INS/GPS and the LGU/GPS at any time t is as follows: delta phi_E，ΔΦ_N，ΔΦ_UWhere subscripts E, N, U represent the east, north, and sky components, respectively. Delta phi_E、ΔΦ_NAnd Δ Φ_UCalculated from equation (2):

$\left\{\begin{array}{c}{\mathrm{\Δ\Φ}}_E=C_{n^{\′}}^n\left(\mathrm{2,3}\right)\\ {\mathrm{\Δ\Φ}}_N=C_{n^{\′}}^n\left(\mathrm{3,1}\right)\\ {\mathrm{\Δ\Φ}}_U=C_{n^{\′}}^n\left(\mathrm{1,2}\right)\end{array}\right.---\left(2\right)$

whereinTo representThe elements of row 2, column 3,to representThe elements of row 3, column 1,to representRow 1, column 2 elements of the matrix.

And 5: and establishing an observation equation and a state equation of the vertical deviation measurement, and extracting the vertical deviation by using a Kalman filtering method with the aid of a global gravity model. The specific implementation method is as follows:

recording the deviation of the true east vertical line of the position of the carrier at the time of measurement as η, the deviation of the true north vertical line as ξ, and recording the deviation values of the east vertical line and the north vertical line calculated by the global gravity model as ηThe EGM2008 global gravity model is selected as the global gravity model, and the calculation program and the use method of the global gravity model can be selected from "http://earth- info.nga.mil/GandG/wgs84/gravitymod/index.htmlThe deviation errors of the east and north vertical lines calculated by the global gravity model of EGM2008 are respectively recorded as η and ξ, and satisfy the relationship of equation (3):

$\left\{\begin{array}{c}\mathrm{\η}=\widehat{\mathrm{\η}}+\mathrm{\δ\η}\\ \mathrm{\ξ}=\widehat{\mathrm{\ξ}}+\mathrm{\δ\ξ}\end{array}\right.---\left(3\right)$

the attitude error angle of INS/GPS output relative to n' system is recorded as v_E，v_N，v_UThe LGU/GPS output has an attitude error angle phi with respect to the n system_E、φ_N、φ_USubscript E, N, U denotes the east, north, and sky components, respectively.

Step 5.1: and establishing a state equation of the vertical deviation measurement.

Selecting a state variable of the plumb line deviation measurement system as phi_E、φ_N、φ_U、η、ξ、_UWherein_URespectively dynamically modeling the state quantities for the equivalent day zero offset of the laser gyro, phi_E、φ_N、φ_UThe following differential equation is satisfied:

$\left\{\begin{array}{c}{\stackrel{\·}{\mathrm{\φ}}}_E=-{\mathrm{\ω}}_{\mathrm{ie}}\cos L\·{\mathrm{\φ}}_U+{\mathrm{\ω}}_{\mathrm{ie}}\sin L\·{\mathrm{\φ}}_N\\ {\stackrel{\·}{\mathrm{\φ}}}_N=-{\mathrm{\ω}}_{\mathrm{ie}}\sin L\·{\mathrm{\φ}}_E\\ {\stackrel{\·}{\mathrm{\φ}}}_U={\mathrm{\ω}}_{\mathrm{ie}}\cos L\·{\mathrm{\φ}}_E-{\mathrm{\ϵ}}_U\end{array}\right.---\left(4\right)$

wherein, ω is_ieIs the rotational angular velocity of the earth, and L is the geographical latitude of the measuring point.

_UModeling is a random constant model, and then:

${\stackrel{\·}{\mathrm{\ϵ}}}_U=0---\left(5\right)$

the statistical models of the deviation errors eta, xi of the east and north plumb lines of the EGM2008 satisfy the following differential equations, respectively:

$\left[\begin{array}{c}{\stackrel{\·}{x}}_E\\ \mathrm{\δ}\stackrel{\·}{\mathrm{\η}}\end{array}\right]=\left[\begin{array}{cc}0& 1\\ -{\mathrm{\ω}}_0^2& -2\mathrm{\ζ}{\mathrm{\ω}}_0\end{array}\right]\left[\begin{array}{c}{x}_E\\ \mathrm{\δ\η}\end{array}\right]+\left[\begin{array}{c}0\\ w_E\end{array}\right]---\left(6\right)$

$\left[\begin{array}{c}{\stackrel{\·}{x}}_N\\ \mathrm{\δ}\stackrel{\·}{\mathrm{\ξ}}\end{array}\right]=\left[\begin{array}{cc}0& 1\\ -{\mathrm{\ω}}_0^2& -2\mathrm{\ζ}{\mathrm{\ω}}_0\end{array}\right]\left[\begin{array}{c}{x}_N\\ \mathrm{\δ\ξ}\end{array}\right]+\left[\begin{array}{c}0\\ w_N\end{array}\right]---\left(7\right)$

wherein x_E、x_NIs an intermediate state quantity; omega₀Is the natural frequency, ω, of the statistical model₀With a fixed relationship omega to the speed of movement V of the carrier₀Zeta is damping coefficient of the statistical model and can be set to be 0.05, w_EAnd w_NProcess noise, w, of the statistical model for the east and north plumb deviations, respectively_E～(0,q_E)，w_N～(0,q_N) I.e. w_EAnd w_NRespectively obey a Gaussian distribution with a mean value of 0 and a variance of q_EAnd q is_NIts typical value may be set to q_E＝(3″)²×(4ζω₀)，q_N＝(3″)²×(4ζω₀)。

The state space vector is selected as:

x＝[φ_E,φ_N,φ_U,_U,x_E,x_N,η,ξ]^T(8)

writing equations (4) - (7) as a unified equation of state is of the form:

$\stackrel{\·}{x}=F\·x+w---\left(9\right)$

wherein, the matrix $F=\left[\begin{array}{cccccccc}0& {\mathrm{\ω}}_{\mathrm{ie}}\sin L& -{\mathrm{\ω}}_{\mathrm{ie}}\cos L& 0& 0& 0& 0& 0\\ -{\mathrm{\ω}}_{\mathrm{ie}}\sin L& 0& 0& 0& 0& 0& 0& 0\\ {\mathrm{\ω}}_{\mathrm{ie}}\cos L& 0& 0& -1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 1\\ 0& 0& 0& 0& -{\mathrm{\ω}}_0^2& 0& -2\mathrm{\ζ}{\mathrm{\ω}}_0& 0\\ 0& 0& 0& 0& 0& -{\mathrm{\ω}}_0^2& 0& -2\mathrm{\ζ}{\mathrm{\ω}}_0\end{array}\right],$

w＝[0,0,0,0,0,0,w_E,w_N]^T. Equation (9) is the state equation of the perpendicular deviation measurement system. The process noise w — (0, Q), i.e., w follows a gaussian distribution with a mean of 0 and a variance of Q. $Q=\left[\begin{array}{cc}{O}_{6\×6}& O_{6\×2}\\ O_{2\×6}& \mathrm{diag}(q_E,q_N)\end{array}\right],$ Wherein O is_6×6、O_2×6、O_6×2A zero matrix, diag (q) representing 6 rows and 6 columns, 2 rows and 6 columns, 6 rows and 2 columns, respectively_E,q_N) Is represented by q_EAnd q is_NIs a diagonal matrix of diagonal elements.

Discretizing the state equation to obtain:

x_i＝M_i/i-1·x_i-1+w_i(10)

x_idenotes x is at t_iSampled value of time, M_i/i-1Is t_i-1To t_iState transition matrix of time, and M_i/i-1≈I_8×8+Δt·F_iIn which I_8×8Is an 8-dimensional identity matrix, F_iFor the matrix F at t_iThe sampled value at the time, Δ t, is the sampling interval. w is a_iObedience mean 0, variance Q_iGaussian distribution of (Q)_i＝Δt²·Q。

Step 5.2: and establishing an observation equation of the vertical deviation measurement.

And establishing the following observation equation by taking the difference between the attitude angles output by the INS/GPS and the LGU/GPS as an observed quantity:

$\left\{\begin{array}{c}{\mathrm{\Δ\Φ}}_E=\mathrm{\ξ}-{\mathrm{\φ}}_E+v_E\\ {\mathrm{\Δ\Φ}}_N=-\mathrm{\η}-{\mathrm{\φ}}_N+v_N\end{array}\right.---\left(11\right)$

substituting the formula (3) for the formula (11) to obtain the following product:

$\left\{\begin{array}{c}{\mathrm{\Δ\Φ}}_E-\widehat{\mathrm{\ξ}}=\mathrm{\δ\ξ}-{\mathrm{\φ}}_E+v_E\\ {\mathrm{\Δ\Φ}}_N+\widehat{\mathrm{\η}}=-\mathrm{\δ\η}-{\mathrm{\φ}}_N+v_N\end{array}\right.---\left(12\right)$

selecting a new observation vector as follows:

$y={[{\mathrm{\Δ\Φ}}_E-\widehat{\mathrm{\ξ}},{\mathrm{\Δ\Φ}}_N+\widehat{\mathrm{\η}}]}^T---\left(13\right)$

rewrite equation (12) to:

y＝H·x+v (14)

wherein the matrix $H=\left[\begin{array}{cccccccc}-1& 0& 0& 0& 0& 0& 0& 1\\ 0& -1& 0& 0& 0& 0& -1& 0\end{array}\right],$ v＝[v_E,v_N]^T。

Equation (14) is an observation equation of the perpendicular deviation measurement system. INS/GPS attitude error angle v_E、v_NObserved noise, v, of the east and north observed components, respectively_E～(0,r_E)，v_N～(0,r_N) I.e. v_EAnd v_NRespectively obey a Gaussian distribution with a mean value of 0 and a variance of r_EAnd r_N. The observation noise variance matrix is based on the observation noise vector v (0, R), i.e., v obeys a Gaussian distribution with mean 0 and variance R $r=\left[\begin{array}{cc}{r}_E& 0\\ 0& r_N\end{array}\right].$ r_EAnd r_NMay be set to 0.5 ".

Discretizing the observation equation to obtain:

y_i＝H_i·x_i+v_i(15)

y_i、v_i、H_iy, v, H at t_iSampled value of time, v_iObedience mean 0, variance R_iOf Gaussian distribution of R_i＝Δt²·R，H_i＝H。

Step 5.3: and (3) estimating the state vector by using a Kalman filtering algorithm according to the state equation and the observation equation of the discrete form of the perpendicular deviation measurement system described by the equations (10) and (15).

First, setting the initial value of the state vectorAnd estimating a noise covariance matrix P₀The initial value of (c).Can be set as an 8-dimensional zero vector; p₀Set as a diagonal matrix, can be set as a typical value:

P₀＝diag((2″)²，(2″)²，(5″)²，(0.003″)²，(1″)²，(1″)²，(1″)²，(1″)²)。

and then updating the state vector and the covariance matrix estimated value by a Kalman filtering iterative algorithm in the formula (16).

$\left\{\begin{array}{c}{\hat{x}}_{i/i-1}=M_{i/i-1}\·{\hat{x}}_{i-1}\\ P_{i/i-1}=M_{i/i-1}\·P_{i-1}{M}_{i/i-1}^T+Q_i\\ K_i=P_{i/i-1}\·H_i^T\·{[H_i\·P_{i/i-1}\·H_i^T+R_i]}^{-1}\\ P_i=(I-K_i\·H_i)\·P_{i/i-1}\·{[I-K_i\·H_i]}^{-1}+K_i\·P_{i/i-1}\·K_i^T\\ {\hat{x}}_i={\hat{x}}_{i/i-1}+K_i\·(y_i-H_i\·{\hat{x}}_{i/i-1})\end{array}\right.---\left(16\right)$

WhereinIs x_iIs determined by the estimated value of (c),is x_iOne step of predicting value of, P_iIs x_iIs estimated by the covariance matrix, P_i/i-1Is P_iOne step of predicting value of, K_iFor filter gain]^-1The symbol represents the inversion operation of the matrix.

t_iThe estimates of the errors η, ξ for the east and north plumb line deviations at the times areFrom t_iState vector estimation at timeThe following are given:

$\left\{\begin{array}{c}\mathrm{\δ}{\widehat{\mathrm{\η}}}_i={\hat{x}}_i\left(7\right)\\ \mathrm{\δ}{\widehat{\mathrm{\ξ}}}_i={\hat{x}}_i\left(8\right)\end{array}\right.---\left(17\right)$

step 5.4: calculating a vertical deviation measurement value:

t_itime east and north plumb deviation measurement η_i、ξ_iCalculated from equations (18) and (19), respectively:

${\mathrm{\η}}_i={\widehat{\mathrm{\η}}}_i+\mathrm{\δ}{\widehat{\mathrm{\η}}}_i---\left(18\right)$

${\mathrm{\ξ}}_i={\widehat{\mathrm{\ξ}}}_i+\mathrm{\δ}{\widehat{\mathrm{\ξ}}}_i---\left(19\right)$

whereinAnd calculating the deviation values of the east-oriented vertical line and the north-oriented vertical line respectively obtained by the global gravity model.

Claims (4)

1. A vertical line deviation dynamic measurement method based on combination of an INS and a GPS is implemented based on a device which comprises an inertial navigation system (1), a GPS antenna (2), a GPS receiver (3) and a data processing computer (4), wherein the inertial navigation system (1), the GPS antenna (2) and the GPS receiver (3) form an INS/GPS attitude measurement subsystem, the inertial navigation system (1) comprises three orthogonally-installed laser gyros which are called laser gyro assemblies, the GPS antenna (2) and the GPS receiver (3) form an LGU/GPS attitude measurement subsystem, and the three laser gyros are all communicated with the GPS receiver (3); the inertial navigation system (1), the GPS antenna (2) and the GPS receiver (3) are fixedly connected and installed on a measurement carrier (5), the GPS antenna (2) is communicated with the GPS receiver (3), the inertial navigation system (1) and the GPS receiver (3) are connected with a data processing computer (4) through data lines, measurement data of the inertial navigation system (1) and the GPS receiver (3) are transmitted into the data processing computer (4) through the data lines, and calculation of vertical line deviation is completed in the data processing computer (4); the method is characterized by comprising the following steps:

1) an INS/GPS attitude measurement subsystem is constructed, and alignment is carried out for more than 8 hours:

starting the INS/GPS attitude measurement subsystem to carry out alignment for more than 8 hours, realizing combined attitude measurement by the INS/GPS attitude measurement subsystem in the alignment process, and continuously outputting an attitude matrix of the INS coordinate system relative to the computed navigation coordinate system by the INS/GPS attitude measurement subsystem in the whole measurement processNo vertical deviation data is output in the alignment process;

2) after the INS/GPS alignment is finished, starting an LGU/GPS attitude measurement subsystem, and initializing an attitude matrix for the LGU/GPS by using the attitude matrix output after the INS/GPS alignment in the first step:

note t_iI is more than or equal to 1 and less than or equal to N at the moment corresponding to the ith measurement sampling point, N is the total number of sampling points in the whole measurement process, t_iThe attitude matrix output by the INS/GPS at the moment isLGU/GPS output attitude ofThe superscript n represents the true navigation coordinate system, and the moment when the LGU/GPS attitude measurement subsystem is started is recorded as t₀The initialization method of LGU/GPS attitude matrix is order

3) And the LGU/GPS attitude measurement subsystem carries out attitude updating:

the LGU/GPS attitude measurement subsystem carries out attitude updating, and meanwhile, the data processing computer stores attitude matrixes output by INS/GPS at all measurement sampling momentsAnd LGU/GPS output attitude matrix

4) Using t in the whole measurement process_iAttitude matrix output by time INS/GPSAnd LGU/GPS output attitude matrixCalculating coordinate rotation matricesAnd calculates the difference between the corresponding attitude angles of its outputs:

the calculation method comprises the following steps:

$C_{n^{\′}\left(i\right)}^{n\left(i\right)}=C_{b\left(i\right)}^{n\left(i\right)}\·{\[C_{b\left(i\right)}^{n^{\′}\left(i\right)}\]}^T---\left(1\right)$

wherein]^TRepresents transpose of matrix, sign-represents matrix multiplication;

simplifying the attitude matrix from n' system to n system obtained by calculating at any time t asThe difference between the three attitude angles output by the INS/GPS and the LGU/GPS at any time t is as follows: delta phi_E，ΔΦ_N，ΔΦ_UWherein the subscript E, N, U denotes the east, north and sky components, respectively, Δ Φ_E、ΔΦ_NAnd Δ Φ_UCalculated from the following formula:

$\left\{\begin{array}{c}{\mathrm{\Δ\Φ}}_E=C_{n^{\′}}^n(2,3)\\ {\mathrm{\Δ\Φ}}_N=C_{n^{\′}}^n(3,1)\\ {\mathrm{\Δ\Φ}}_U=C_{n^{\′}}^n(1,2)\end{array}\right.---\left(2\right)$

whereinTo representThe elements of row 2, column 3,to representThe elements of row 3, column 1,to representRow 1, column 2 elements of the matrix;

5) the difference delta phi of the three attitude angles output by the INS/GPS and the LGU/GPS_E，ΔΦ_N，ΔΦ_UFor the observed quantity, the vertical deviation is extracted by a Kalman filtering method with the assistance of a global gravity model by establishing an observation equation and a state equation of the vertical deviation measurement:

the deviation of the true east perpendicular line of the position of the measuring carrier at the moment of t is recorded as η, the deviation of the true north perpendicular line is recorded as ξ, and the deviation values of the east perpendicular line and the north perpendicular line obtained by calculation of the global gravity model are respectively recorded asThe deviation errors of the east and north vertical lines calculated by the global gravity model are respectively recorded as η and ξ, and have the following relationship:

$\left\{\begin{array}{c}\mathrm{\η}=\widehat{\mathrm{\η}}+\mathrm{\δ}\mathrm{\η}\\ \mathrm{\ξ}=\widehat{\mathrm{\ξ}}+\mathrm{\δ}\mathrm{\ξ}\end{array}\right.---\left(3\right)$

the attitude error angle of INS/GPS output relative to n' system is recorded as v_E，v_N，v_UThe LGU/GPS output has an attitude error angle phi with respect to the n system_E、φ_N、φ_U(ii) a The specific implementation method comprises the following steps:

5.1) establishing a state equation of the perpendicular deviation measurement:

selecting a state variable of the plumb line deviation measurement system as phi_E、φ_N、φ_U、η、ξ、_UWherein_URespectively carrying out dynamic modeling on the state variables for the equivalent day zero offset of the laser gyro, phi_E、φ_N、φ_UThe following differential equation is satisfied:

$\left\{\begin{array}{c}{\stackrel{\·}{\mathrm{\φ}}}_E=-{\mathrm{\ω}}_{ie}\cos L\·{\mathrm{\φ}}_U+{\mathrm{\ω}}_{ie}\sin L\·{\mathrm{\φ}}_N\\ {\stackrel{\·}{\mathrm{\φ}}}_N=-{\mathrm{\ω}}_{ie}\sin L\·{\mathrm{\φ}}_E\\ {\stackrel{\·}{\mathrm{\φ}}}_U={\mathrm{\ω}}_{ie}\cos L\·{\mathrm{\φ}}_E-{\mathrm{\ϵ}}_U\end{array}\right.---\left(4\right)$

wherein, ω is_ieThe rotational angular velocity of the earth is shown, and L is the geographical latitude of the measuring point;

_Umodeling is a random constant model, and then:

${\stackrel{\·}{\mathrm{\ϵ}}}_U=0---\left(5\right)$

the statistical models of the deviation errors eta and xi of the east and north vertical lines of the global gravity model are respectively given by the following two formulas:

wherein x_E、x_NIs an intermediate state variable, ω₀Is the natural frequency, ω, of the statistical model₀With a fixed relationship omega to the speed of movement V of the measurement carrier₀＝2π×V/10000，Is the damping coefficient, w, of the statistical model_EAnd w_NProcess noise, w, of the statistical model for the east and north plumb deviations, respectively_E～(0,q_E)，w_N～(0,q_N) I.e. w_EAnd w_NRespectively obey a Gaussian distribution with a mean value of 0 and a variance of q_EAnd q is_N；

The state space vector of the selected vertical deviation measurement system is as follows:

x＝[φ_E,φ_N,φ_U,_U,x_E,x_N,η,ξ]^T(8)

writing equations (4) - (7) as a unified equation of state is of the form:

$\stackrel{\·}{x}=F\·x+w---\left(9\right)$

wherein, the matrixw＝[0,0,0,0,0,0,w_E,w_N]^TEquation (9) is the state equation of the vertical deviation measurement system, the process noise w-0, Q, i.e. the gaussian distribution with the mean value of 0 and the variance of Q,wherein O is_6×6、O_2×6、O_6×2A zero matrix representing 6 rows and 6 columns, 2 rows and 6 columns, 6 rows and 2 columns, respectively, diag [ q ]_E,q_N]Is represented by q_EAnd q is_NA diagonal matrix being diagonal elements;

discretizing equation of state (9) yields:

x_i＝M_i/i-1·x_i-1+w_i(10)

x_idenotes x is at t_iSampled value of time, M_i/i-1Is t_i-1To t_iState transition matrix of time, M_i/i-1≈I_8×8+Δt·F_iIn which I_8×8Is an 8-dimensional identity matrix, F_iFor the matrix F at t_iSample value at time, Δ t is the sampling interval, w_iObedience mean 0, variance Q_iGaussian distribution of (Q)_i＝Δt²·Q；

5.2) establishing an observation equation of the vertical deviation measurement:

the difference delta phi between the attitude angles output by INS/GPS and LGU/GPS_E，ΔΦ_NFor the observed quantity, the following observation equation is established:

$\left\{\begin{array}{c}{\mathrm{\Δ\Φ}}_E=\mathrm{\ξ}-{\mathrm{\φ}}_E+v_E\\ {\mathrm{\Δ\Φ}}_N=-\mathrm{\η}-{\mathrm{\φ}}_N+v_N\end{array}\right.---\left(11\right)$

substituting formula (3) into formula (11) to obtain:

$\left\{\begin{array}{c}{\mathrm{\Δ\Φ}}_E-\widehat{\mathrm{\ξ}}=\mathrm{\δ}\mathrm{\ξ}-{\mathrm{\φ}}_E+v_E\\ {\mathrm{\Δ\Φ}}_N+\widehat{\mathrm{\η}}=-\mathrm{\δ}\mathrm{\η}-{\mathrm{\φ}}_N+v_N\end{array}\right.---\left(12\right)$

selecting a new observation vector as follows:

$y={\[{\mathrm{\Δ\Φ}}_E-\widehat{\mathrm{\ξ}},{\mathrm{\Δ\Φ}}_N+\widehat{\mathrm{\η}}\]}^T---\left(13\right)$

rewrite equation (12) to:

y＝H·x+v (14)

wherein the matrixv＝[v_e，v_N]^T；

Equation (14) is the observation equation of the vertical deviation measurement system, and the INS/GPS attitude error angle v_E、v_NObserved noise, v, of the east and north observed components, respectively_E～(0,r_E)，v_N～(0,r_N) I.e. v_EAnd v_NRespectively obey a Gaussian distribution with a mean value of 0 and a variance of r_EAnd r_NThe observation noise variance matrix is the Gaussian distribution of the observation noise v (0, R), i.e., v obeys a mean of 0 and a variance of R

Discretizing the observation equation to obtain:

y_i＝H_i·x_i+v_i(15)

y_i、v_i、H_iy, v, H at t_iSampled value of time, v_iObedience mean 0, variance R_iOf Gaussian distribution of R_i＝Δt²·R，H_i＝H；

5.3) estimating the state vector by using a Kalman filtering algorithm according to the state equation and the observation equation of the discrete form of the perpendicular deviation measurement system described by the equations (10) and (15):

the iterative algorithm of the Kalman filter algorithm is as follows:

$\left\{\begin{array}{c}{\hat{x}}_{i/i-1}=M_{i/i-1}\·{\hat{x}}_{i-1}\\ P_{i/i-1}=M_{i/i-1}\·P_{i-1}\·M_{i/i-1}^T+Q_i\\ K_i=P_{i/i-1}\·H_i^T\·{\[H_i\·P_{i/i-1}\·H_i^T+R_i\]}^{-1}\\ P_i=(I-K_i\·H_i)\·P_{i/i-1}\·{\[I-K_i\·H_i\]}^{-1}+K_i\·P_{i/i-1}\·K_i^T\\ {\hat{x}}_i={\hat{x}}_{i/i-1}+K_i\·(y_i-H_i\·{\hat{x}}_{i/i-1})\end{array}\right.---\left(16\right)$

whereinIs x_iIs determined by the estimated value of (c),is x_iOne step of predicting value of, P_iIs x_iIs estimated by the covariance matrix, P_i/i-1Is P_iOne step of predicting value of, K_iFor filter gain]^-1Sign represents the inversion operation of the matrix;

t_ithe estimates of the errors η, ξ for the east and north plumb line deviations at the times areFrom t_iState vector estimation at timeThe following are given:

$\left\{\begin{array}{c}\mathrm{\δ}{\widehat{\mathrm{\η}}}_i={\hat{x}}_i\left(7\right)\\ \mathrm{\δ}{\widehat{\mathrm{\ξ}}}_i={\hat{x}}_i\left(8\right)\end{array}\right.---\left(17\right)$

wherein,respectively representing state vector estimation valuesThe 7 th and 8 th elements of (1);

5.4) calculating the vertical deviation measurement value:

t_itime east and north plumb deviation measurement η_i、ξ_iCalculated from equations (18) and (19), respectively:

${\mathrm{\η}}_i={\widehat{\mathrm{\η}}}_i+\mathrm{\δ}{\widehat{\mathrm{\η}}}_i---\left(18\right)$

${\mathrm{\ξ}}_i={\widehat{\mathrm{\ξ}}}_i+\mathrm{\δ}{\widehat{\mathrm{\ξ}}}_i---\left(19\right)$

whereinT respectively calculated for global gravity model_iThe time measures the value of the east and north vertical deviation at the point.

2. A method for dynamically measuring vertical deviation based on a combination of INS and GPS as claimed in claim 1, wherein: the measuring carrier (5) is a measuring ship or a measuring vehicle.

3. A method for dynamically measuring vertical deviation based on a combination of INS and GPS as claimed in claim 1, wherein: and the INS adopts a single-shaft rotary laser gyro inertial navigation system.

4. A method for dynamically measuring vertical deviation based on a combination of INS and GPS as claimed in claim 1, wherein: the global gravity model adopts an EGM2008 global gravity model.