Low-orbit satellite strapdown inertial navigation / radar altimeter combined navigation method and terminal device

Abstract

The application discloses a low-orbit satellite strapdown inertial navigation / radar altimeter combined navigation method and a terminal device. Firstly, the non-spherical perturbation J4 term of the earth and atmospheric resistance are considered to accurately establish an orbit dynamics model, radar altimeter measurement data is introduced to establish a radar altimeter measurement equation, atmospheric resistance is measured by a satellite-borne accelerometer, and an accelerometer measurement equation is established accordingly, and an extended Kalman filter is used to estimate position and speed. The application is a navigation method combining a radar altimeter and a strapdown inertial navigation, can effectively reduce navigation error, and improve autonomous navigation capability under a low-orbit satellite dynamics environment.

Description

Technical Field

[0001] This invention relates to the field of integrated navigation including inertial navigation, and more specifically to a strapdown inertial navigation / radar altimeter integrated navigation method and terminal equipment for low-Earth orbit satellites. Background Technology

[0002] Low Earth orbit (LEO) satellites generally refer to artificial Earth satellites with an orbital period of 84-127 minutes and an orbital altitude of less than 2000 km. Due to their low orbital altitude, they can perform high-precision geodesy, resource surveys, and meteorological observations by collecting Earth observation data. Accurately obtaining satellite positions is a prerequisite for scientific research using LEO satellite data. Currently, the commonly used navigation methods are ground-based radio navigation and satellite navigation. Ground-based radio navigation has many shortcomings, such as long latency, low accuracy, and poor anti-interference performance, and its dependence on ground stations also limits its further development. When using GPS or BeiDou satellite navigation, there are shortcomings such as the signal being easily interfered with, which limits the application scenarios.

[0003] Inertial navigation, as a classic fully autonomous navigation system, boasts numerous advantages such as high short-term accuracy, continuous output, strong anti-interference capabilities, and comprehensive navigation information, leading to its widespread application. However, even expensive, high-precision inertial navigation equipment inevitably suffers from gyroscope and accelerometer drift, causing navigation errors to accumulate over time. For low-Earth orbit (LEO) satellites, due to their lower orbital altitude, they still experience slight atmospheric drag during flight, resulting in deviations between actual accelerometer measurements and ideal values. This navigation error divergence is more pronounced compared to medium- and high-Earth orbit (MEO) satellites. This gradually increasing error over time is unacceptable for LEO satellites requiring long-term operation, necessitating the use of other navigation methods to correct inertial navigation errors. Radar altimeters, as active distance sensors, offer advantages such as a large altitude measurement range, high accuracy, low power consumption, and small size, providing precise altitude measurement information for flight platforms. Liu Jun et al. proposed an autonomous navigation method for radar altimeters based on UKF (Liu Jun, Han Chao. Autonomous orbit determination of radar altimeters based on UKF [J]. Journal of Beijing University of Aeronautics and Astronautics, 2006(08):889-893.); Videmsek et al. (Videmsek A, Haag M, Bleakley T. Radar Altimeter Aiding of GNSS for Precision Approach and Landing of RPA [C] / / 2019 Integrated Communications, Navigation and Surveillance Conference (ICNS). 2019. and Videmsek A, Haag M. On-Board GPS Augmentation through RADAR Altimeter Aiding for Precision Approach and Landing of UAS [C] / / 2020 European Navigation The conference (ENC).2020) proposed a radar altimeter-assisted GNSS navigation method, which combines radar altimeter measurement data with a terrain elevation database to obtain additional distance measurement information to assist GNSS navigation. The above-mentioned literature has studied the radar altimeter navigation method, but it has given little consideration to factors such as the non-spherical perturbation of the earth and atmospheric drag errors in the modeling process. The dynamic model is relatively simple and not suitable for low-Earth orbit satellites. Summary of the Invention

[0004] The technical problem to be solved by the present invention is to provide a low-Earth orbit satellite strapdown inertial navigation / radar altimeter integrated navigation method and terminal equipment to improve the accuracy of low-Earth orbit satellite autonomous navigation, in order to address the shortcomings of the existing technology.

[0005] To solve the above-mentioned technical problems, the technical solution adopted by the present invention is: a low-orbit satellite strapdown inertial navigation / radar altimeter integrated navigation method, comprising the following steps:

[0006] 1) Given the initial state information of the low-Earth orbit satellite, let the J2000.0 geocentric protocol inertial coordinate system be the i-frame and the low-Earth orbit satellite body coordinate system be the b-frame. The initial state information includes the initial time t0 and the terminal time t0. end The directional cosine matrix C of the i-system and the b-system i b The initial position r0, initial velocity v0, initial position error δr0, initial velocity error δv0, and zero offset error b of the low-orbit satellite in the i-frame. f ;

[0007] Let counter k = 0, t k The state of the low-orbit satellite in the i-system at that moment Error status =0 (0 vector), given the initial state error root mean square deviation P k Initial noise variance matrix Q k Given the measurement noise matrix R k r0 and v0 are both 3-dimensional column vectors; X k It is a 6-dimensional column vector;

[0008] 2) Based on the dynamic model considering the non-spherical perturbation J4 term of the Earth, according to t k State X at time X k Predict the motion state at the next moment to obtain t k+1 State estimate X at time 1 k+1,k :

[0009] X k+1,k =Φ k+1,k X k +Γ k b f ;

[0010] Where, Φ k+1,k For t k Time to t k+1 The state transition matrix at time Φ k+1,k =I 6×6 +A·Δt,Δt=t k+1 -t k , I is the identity matrix, O is the zero matrix; Γ k For tk The noise driving matrix at time Γ k =0.5(I 6×6 +Φ k+1,k )·B·Δt, for Transpose of; μ is the gravitational constant, and x, y, z are the gravitational constants. k The three components, r is r k The magnitudes of J2, J3, and J4 are the 2nd, 3rd, and 4th order zonal harmonic terms of the Earth's nonspherical perturbation, respectively, and R... e It is the semi-major axis of the Earth's ellipsoid; c2, c 31 c 33 c 41 It is a diagonal matrix;

[0011] 3) According to t k Error state δX at time 1 k and the mean square error P k The error state and mean square error of the next time step are predicted to obtain t. k+1 Error state estimate δX at time 1 k+1,k and the mean square error estimate P k+1,k The prediction formula is:

[0012] δX k+1,k =Φ k+1,k δX k +Γ k b f ;

[0013]

[0014] in, For Φ k+1,k transpose, For Γ k Transpose of;

[0015] 4) According to t k+1 Measurement equations are established using radar altimeter and accelerometer data at each time point, and the error state δX at the current time point is estimated. k+1,k and the mean square error estimate P k+1,k Perform a measurement update to obtain the updated error state δX. k+1 and the mean square error P k+1 :

[0016] δX k+1 =δX k+1,k +K k+1 (Z k+1 -H k+1δX k+1,k );

[0017]

[0018] Among them, Z k+1 For t k+1 Measurement of time, H k+1 For t k+1 Measurement matrix at time, R k+1 For t k+1 Measurement noise matrix at time step For t k+1 Filter gain at time 10:00 For K k+1 Transpose of;

[0019] 5) Based on the updated error state δX k+1 For t k+1 State estimate X at time 1 k+1 ,k is compensated to obtain t k+1 The state of a low-Earth orbit satellite in the i-system at any given time. k+1 :

[0020] X k+1 =X k+1,k -δX k+1 ;

[0021] 6) Determine t k+1 <t end Is it valid?

[0022] If so, then the X state of the low-orbit satellite in the i-system will be... k+1 Assigned to X k δX k =0, updated mean squared error P k+1 Assign to P k , t k+1 Measurement noise matrix R at time step k+1 Assign to R k , t k+1 Measurement matrix H at time k+1 Assigned to H k Increment the value of k by 1 and return to step 2);

[0023] If not, navigation ends, output the terminal's state X and the root mean square error P at that time:

[0024] X = X k+1 ;

[0025] P = P k+1 .

[0026] This invention considers the J4 term of Earth's non-spherical perturbation and the influence of atmospheric drag, establishing a more accurate orbital dynamics model. A geoid model is used to replace the actual Earth surface, and the difference between the altitude measured by a radar altimeter and the calculated altitude obtained from inertial navigation is used as a measurement to derive the radar altimeter measurement equation. An accelerometer measurement equation is established based on atmospheric drag measured by an accelerometer, and an extended Kalman filter is used to estimate the position and velocity of low-Earth orbit (LEO) satellites. The dynamics model is well-considered, significantly improving the accuracy of autonomous navigation for LEO satellites.

[0027] In this invention, in:

[0028]

[0029] c2 = diag(1.5, 1.5, 4.5), c 31 =diag(7.5,7.5,15), c 33 =diag(0,0,1.5), c 41 =diag(1.875,1.875,9.375), c 42 =diag(26.25,26.25,43.75).

[0030] In this invention, t k+1 The measurement of time is:

[0031] Z k+1 =H k+1 X k+1,k +υ;

[0032] Where υ is the measurement error, υ = [υ h υ a ] T υ h For radar altimeter measurement error, υ a To account for the measurement error of the accelerometer, Z h,k+1 For t k+1 Measurement of the radar altimeter equation, Z a,k+1 For t k+1 The measurement of the equation of acceleration at any given time, Z k+1 =[Z h,k+1 Z a,k+1 ] T .

[0033] The measurement method of this invention takes into account the measurement errors of radar altimeters and accelerometers, thereby improving the accuracy of measurement calculations and further enhancing the accuracy of low-orbit satellite navigation.

[0034] In this invention, t k+1Measurement of radar altimetry equation Where h h,k+1 For t k+1 Elevation data measured by radar altimeter at all times. For t k+1 Time based on position estimate r k+1,k Speed ​​estimate v k+1,k The calculated elevation data. In this invention, the equations on both sides of the radar altimeter use a geoid model instead of the actual Earth surface when calculating elevation, which reduces the impact of measurement errors compared to directly using an ellipsoidal model.

[0035] In this invention, the calculated elevation data According to t k+1 Time position estimate r k+1,k The calculation is as follows:

[0036]

[0037] Among them, lr k+1,k For r k+1,k The magnitude of ξ is the compensation value for terrain undulation, lR k+1,k The radius of the reference ellipsoid is the position of the nadir point of the low-Earth orbit satellite. R is the average radius of the Earth, x k+1 y k+1 z k+1 For r k+1,k The three components, n = 0.00344978650678.

[0038] This invention is based on t k+1 Time position estimate r k+1,k The calculated elevation data Compensation for terrain undulations improves navigation accuracy.

[0039] In this invention, t k+1 Measurement of the equation of acceleration measurement at any time Where a k+1 For t k+1 Atmospheric drag experienced by low-Earth orbit satellites at all times. In this invention, the accelerometer measures atmospheric drag, and based on this, an accelerometer measurement equation is established, reducing the impact of atmospheric drag on navigation accuracy.

[0040] In this invention, a k+1 According to t k+1 Estimated velocity v at time step k+1,k The calculation is as follows:

[0041]

[0042] Among them, C D Here, A is the drag coefficient, m is the frontal area of ​​the low-Earth orbit satellite, ρ is the atmospheric density at the location of the low-Earth orbit satellite, and lv is the atmospheric density at the location of the low-Earth orbit satellite. k+1,k For v k+1,k Size.

[0043] As an inventive concept, the present invention also provides a terminal device comprising a processor and a memory; the memory storing computer programs / instructions; the processor executing the computer programs / instructions stored in the memory; the computer programs / instructions being configured to implement the steps of the method of the present invention.

[0044] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0045] 1. This invention takes into account the J4 term of the Earth's non-spherical perturbation and the influence of atmospheric drag, and establishes a more accurate orbital dynamics model, reducing the impact of model errors on navigation accuracy;

[0046] 2. This invention introduces a radar altimeter measurement equation, which uses a geoid model instead of the actual Earth surface when calculating elevation, thus reducing the impact of measurement errors compared to directly using an ellipsoid model.

[0047] 3. Spaceborne accelerometers operate in a weightless state in space. For the research object of this invention, the accelerometer measures atmospheric drag, and based on this, an accelerometer measurement equation was established to reduce the impact of atmospheric drag on navigation accuracy. Attached Figure Description

[0048] Figure 1 This is a schematic diagram of radar altimeter measurement.

[0049] Figure 2 This is a schematic diagram of the strapdown inertial navigation / radar altimeter integrated navigation system. Detailed Implementation

[0050] This invention addresses the problem of high-precision autonomous navigation for low-Earth orbit (LEO) satellites by proposing a strapdown inertial navigation system (INS) / radar altimeter integrated navigation method. This method utilizes radar altimeter measurement data to correct the INS error, suppressing INS error divergence. The basic idea is as follows: First, considering the J4 term of Earth's non-spherical perturbation and the influence of atmospheric drag, an accurate orbital dynamics model is established. Second, based on the radar altimeter measurement principle, a geoid model is used to replace the actual Earth surface, deriving the elevation calculation formula and establishing the radar altimeter measurement equation. Finally, the extended Kalman filter method is used to estimate and correct the position and velocity errors of the LEO satellite.

[0051] This invention mainly includes the following steps:

[0052] Step 1: Provide the initial state information of the low-Earth orbit satellite. Let the J2000.0 geocentric protocol inertial coordinate system be the i-frame, and the low-Earth orbit satellite body coordinate system be the b-frame. The state information includes the initial time t0 and the terminal time t... end The directional cosine matrix C of the i-system and the b-system i b The initial position r0, initial velocity v0, initial position error δr0, initial velocity error δv0, and zero offset error b of the low-orbit satellite in the i-frame. f .

[0053] Step 2: Set counter k = 0, t k The state of the low-orbit satellite in the i-system at that moment Error status Given the initial state error mean squared P k Initial noise variance matrix Q k Given the measurement noise matrix R k .

[0054] Step 3: Error state prediction.

[0055] 1) Based on the dynamic model considering the J4 term of the Earth's non-spherical perturbation, according to t k The state variable X at time t is k Predict the motion state at the next moment to obtain t k+1 The time t is obtained k+1 State estimate X at time 1 k+1,k :

[0056] X k+1,k =Φ k+1,k X k +Γ k b f ;

[0057] Where, Φ k+1,k For t k Time to t k+1 The state transition matrix at time Φ k+1,k =I 6×6 +A·Δt,Δt=t k+1 -t k , I is the identity matrix, O is the zero matrix, and the expression for G is given in the appendix; Γ k For t k The noise driving matrix at time Γ k =0.5(I 6×6 +Φ k+1,k )·B·Δt, for transpose, G 11 ~G 33 The expression is as follows:

[0058]

[0059] Where μ is the gravitational constant, c² = diag(1.5, 1.5, 4.5), c 31 =diag(7.5,7.5,15), c 33 =diag(0,0,1.5), c 41 =diag(1.875,1.875,9.375),c 42 =diag(26.25,26.25,43.75), where x, y, and z are r k The three components, r is r k The magnitudes of J2, J3, and J4 are the 2nd, 3rd, and 4th order zonal harmonic terms of the Earth's nonspherical perturbation, respectively, and R... e It is the semi-major axis of the Earth's ellipsoid; in this embodiment of the invention, the value of the diagonal matrix is ​​derived from the calculation formula of the Earth's gravitational acceleration, reference: Zhang Hongbo. Theory and Method of Spacecraft Orbital Mechanics [M]. National Defense Industry Press, 2015.

[0060] 2) According to t k Error state δX at time 1 k and the mean square error P k The error state and mean square error of the next time step are predicted to obtain the result for time step t. k+1 Error state estimate δX at time 1 k+1,k and the mean square error estimate P k+1,k The prediction formula is:

[0061] δX k+1,k =Φ k+1,k δX k +Γ k b f ;

[0062]

[0063] in, For Φ k+1,k transpose, For Γ k The transpose of .

[0064] Step 4: Measurement Update.

[0065] 1) According to t k+1 Time-state estimate X k+1,k The reference ellipsoid radius lR for calculating the nadir position of a low-Earth orbit satellite k+1,k Solving elevation data and atmospheric drag a k+1 :

[0066]

[0067]

[0068]

[0069] Where R is the average radius of the Earth, x k+1 y k+1 z k+1 For r k+1,k The three components, lr k+1,k For r k+1,k The magnitude of ξ is the compensation value for terrain undulation. n = 0.00344978650678; C D Here, A is the drag coefficient, m is the frontal area of ​​the low-Earth orbit satellite, ρ is the atmospheric density at the location of the low-Earth orbit satellite, and lv is the atmospheric density at the location of the low-Earth orbit satellite. k+1,k For v k+1,k Size;

[0070] 2) According to t k+1 Elevation data h measured by radar altimeter at all times h,k+1 Calculate elevation data h h0,k+1 and atmospheric drag a k+1 Establish separate equations for radar altitude measurement and acceleration measurement;

[0071] Radar altimeter measurement equation:

[0072]

[0073] Accelerometer equation:

[0074]

[0075] Then t k+1 The total measurement equation at time t is:

[0076] Z k+1 =H k+1 X k+1,k +υ;

[0077] Among them, Z k+1 =[Z h,k+1 Z a,k+1 ] T υ=[υ h υ a ] T , υ represents the measurement error, υ = [υ hυ a ] T υ h For radar altimeter measurement error, υ a To account for the measurement error of the accelerometer, Z h,k+1 For t k+1 Measurement of the radar altimeter equation, Z a,k+1 For t k+1 The measurement of the equation of acceleration at any given time, Z k+1 =[Z h,k+1 Z a,k+1 ] T . t k+1 Measurement of radar altimetry equation Where h h,k+1 For t k+1 Elevation data measured by radar altimeter at all times. For t k+1 Time based on position estimate r k+1,k Speed ​​estimate v k+1,k The calculated elevation data.

[0078] 3) The error state δX at the current moment k+1,k and the mean square error P k+1,k Perform a measurement update to obtain the updated error state δX. k+1 and the mean square error P k+1 :

[0079] δX k+1 =δX k+1,k +K k+1 (Z k+1 -H k+1 δX k+1,k );

[0080]

[0081] 4) Based on the updated error state δX k+1 For t k+1 State estimate X at time 1 k+1,k Compensation is performed to obtain t k+1 State X at time X k+1 :

[0082] X k+1 =X k+1,k -δX k+1 .

[0083] Step 5: Determine t k+1 <t end Is it valid?

[0084] If so, then the X state of the low-orbit satellite in system i will be... k+1 Assigned to X k δX k =0, updated mean squared error P k+1 Assign to P k , t k+1 Measurement noise matrix R at time step k+1 Assign to R k , t k+1 Measurement matrix H at time k+1 Assigned to H k Increment the value of k by 1 and return to step 3);

[0085] If not, navigation ends, output the terminal's state X and the root mean square error P at that time:

[0086] X = X k+1 ;

[0087] P = P k+1 .

[0088] The feasibility of the present invention will be verified through simulation experiments.

[0089] A satellite orbit with an altitude of 200 km was simulated. The simulation considered the effects of Earth's non-spherical perturbation and atmospheric drag. The gravitational field model selected was EGM2008 of order 100, and the atmospheric model was NRLMSIR-00, using the atmospheric density of a year with high solar activity. The radar altimeter operated at a frequency of 1 Hz, working for 10 seconds every 1000 seconds, with only the inertial navigation system operating during the remaining time periods.

[0090] S1. The initial flight time corresponds to UTC time 00:00:00 on April 2, 2022, the flight time is 7200s, and the initial time t k =0s, terminal time t end =7200s, the reference coordinate system is selected as J2000.0 geocentric protocol inertial coordinate system, the nominal orbital position velocity of the low-orbit satellite in the i-frame is RV, the initial position velocity [r0; v0] = RV0 + [200; 200; 200; 0.5; 0.5; 0.5], add zero offset b. f =20μg;

[0091] S2. Set counter k = 0, the state of the low-orbit satellite in frame i. Error status Initial state error root mean square deviation P k =diag(200) 2 200 2 200 2 0.5 2 0.5 2 0.5 2), initial noise variance matrix Measurement noise matrix

[0092] S3. Based on the dynamic model considering the non-spherical perturbation of the Earth (J4 term), according to t k The state variable X at time t is k Predict the motion state at the next moment to obtain t k+1 State estimate X at time 1 k+1,k :

[0093] X k+1,k =Φ k+1,k X k +Γ k b f ;

[0094] Where, Φ k+1,k For t k Time to t k+1 The state transition matrix at time Φ k+1,k =I 6×6 +A·Δt,Δt=t k+1 -t k , I is the identity matrix, O is the zero matrix; Γ k For t k The noise driving matrix at time Γ k =0.5(I 6×6 +Φ k+1,k )·B·Δt, for The transpose of .

[0095] S4. According to t k Error state δX at time 1 k and the mean square error P k The error state and mean square error of the next time step are predicted to obtain the result for time step t. k+1 Error state estimate δX at time 1 k+1,k and the mean square error estimate P k+1,k The prediction formula is:

[0096] δX k+1,k =Φ k+1,k δX k +Γ k b f ;

[0097]

[0098] S5. According to t k+1 Time-state estimate X k+1,kThe reference ellipsoid radius lR for calculating the nadir position of a low-Earth orbit satellite k+1,k Solving elevation data and atmospheric drag a k+1 Combined with the elevation data h measured by the radar altimeter h,k+1 Establish the measurement equation. Based on the measurement equation, determine the error state δX at the current time. k+1,k and the mean square error P k+1,k Perform a measurement update to obtain the updated error state δX. k+1 and the mean square error P k+1 :

[0099] δX k+1 =δX k+1,k +K k+1 (Z k+1 -H k+1 δX k+1,k );

[0100]

[0101] S6. Based on the updated error state δX k+1 For t k+1 State estimate X at time 1 k+1,k Compensation is performed to obtain t k+1 State X at time X k+1 :

[0102] X k+1 =X k+1,k -δX k+1 ;

[0103] S7. Determine t k+1 <t end Is it valid?

[0104] If so, then the X state of the low-orbit satellite in system i will be... k+1 Assigned to X k δX k =0, updated mean squared error P k+1 Assign to P k , t k+1 Measurement noise matrix R at time step k+1 Assign to R k , t k+1 Measurement matrix H at time k+1 Assigned to H k Increment the value of k by 1 and return S3;

[0105] If not, navigation ends, based on the terminal's current state X. k+1 Calculate the navigation error δX and the root mean square error P:

[0106] δX=Xk+1 -RV k+1 ;

[0107] P = P k+1 .

[0108] The terminal position error was 157.732m, the root mean square error of the position error was 239.097m, the velocity error was 0.173m / s, and the root mean square error of the velocity error was 0.262m / s. When using pure inertial navigation, the root mean square error of the position error was 6377.108m, and the root mean square error of the velocity error was 7.964m / s. When using the integrated navigation method of this embodiment, error divergence is effectively suppressed, and navigation accuracy is significantly improved.

[0109] Table 1 Comparison of accuracy between pure inertial navigation and integrated navigation

[0110]

Claims

1. A low-Earth orbit satellite strapdown inertial navigation / radar altimeter integrated navigation method, characterized in that, Includes the following steps: 1) Given the initial state information of the low-Earth orbit satellite, let the J2000.0 geocentric protocol inertial coordinate system be denoted as... i The coordinate system of the low-orbit satellite body is... b The initial state information includes the initial time. t 0, Terminal Time t end , i System and b Directional cosine array C i b Low-orbit satellites i Initial position in the system r 0 、 initial velocity v 0 and initial position error δr 0. Initial velocity error δv 0 、 Adding zero bias error b f ; Let the counter k = 0, t k Low-Earth orbit satellites at that moment i System state Error status = 0 Given the initial state error mean square error P k 、 Initial noise variance matrix Q k Given a measurement noise matrix R k ; 2) Based on the dynamic model considering the J4 term of the Earth's non-spherical perturbation, according to t k state of time X k Predict the motion state at the next moment and obtain t k+1 State estimate at time 1 X k+1,k : ； in, for t k Time's up t k+1 The state transition matrix at time t, , , , I It is the identity matrix. O It is a zero matrix; for t k The noise driving matrix at time step, , , for Transpose of; , , μ The gravitational constant is the constant of gravity. x , y , z for r k The three components, r for r k Size, J 2. J 3. J 4 represents the 2nd, 3rd, and 4th order zonal harmonic terms of Earth's nonspherical perturbation. R e It is the semi-major axis of the Earth's ellipsoid; , , , It is a diagonal matrix; 3) According to t k Error state at time δX k and the mean square error P k The error state and mean square error at the next time step are predicted to obtain... t k+1 Error state estimate at time 1 δX k+1,k and the mean square error estimate P k+1,k The prediction formula is: ； ； in, for transpose, for Transpose of; 4) According to t k+1 Measurement equations are established using radar altimeter and accelerometer data at each time point to estimate the error state at the current time point. δX k+1,k and the mean square error estimate P k+1,k Perform a measurement update to obtain the updated error state. δX k+1 and the mean square error P k+1 : ； ； in, for t k+1 Measurement of time, for t k+1 Measurement matrix at time, for t k+1 Measurement noise matrix at time step for t k+1 Filter gain at time 10:00 for Transpose of; 5) Based on the updated error status δX k+1 right t k+1 State estimate at time 1 X k+1,k To receive compensation t k+1 Low-Earth Orbit Satellites i The state in the system X k+1 : ； 6) Judgment t k+1 < t end Is it valid? If so, then the low-orbit satellite will be in i The state in the system X k+1 Assign to X k , δX k = 0, updated mean squared error P k+1 Assign to P k , t k+1 Measurement noise matrix at time step R k+1 Assign to R k , t k+1 Measurement matrix at time H k+1 Assign to H k , k Increment the value by 1, then return to step 2); If not, navigation ends, output the terminal's current status. X and the mean square error P : ； ； t k+1 The measurement of time is: ； in, υ To account for measurement error, , υ h Due to radar altimeter measurement error, υ a To account for the measurement error of the accelerometer, , Z h,k+1 for t k+1 Measurement of radar altimetry equations at all times. Z a,k+1 for t k+1 The measurement of the equation for the accelerometer at any given time. ; t k+1 Measurement of radar altimetry equation ,in h h,k+1 for t k+1 Elevation data measured by radar altimeter at all times. for t k+1 Time based on position estimate r k+1,k 、 Speed ​​estimate v k+1,k The calculated elevation data; The calculated elevation data according to t k+1 Time position estimate r k+1,k The calculation is as follows: ； in, lr k+1,k for r k+1,k Size, ξ This is the compensation value for terrain undulations. lR k+1,k The radius of the reference ellipsoid is the position of the nadir point of the low-Earth orbit satellite. ; R The average radius of the Earth x k+1 , y k+1 , z k+1 for r k+1,k The three components, , n =0.00344978650678。 2. The low-orbit satellite strapdown inertial navigation / radar altimeter integrated navigation method according to claim 1, characterized in that, ,in: 。 3. The low-orbit satellite strapdown inertial navigation / radar altimeter integrated navigation method according to claim 1, characterized in that, ， ， ， ， 。 4. The low-orbit satellite strapdown inertial navigation / radar altimeter integrated navigation method according to claim 1, characterized in that, t k+1 Measurement of the equation of acceleration measurement at any time ,in a k+1 for t k+1 Atmospheric drag experienced by low-Earth orbit satellites at all times.

5. The low-Earth orbit satellite strapdown inertial navigation / radar altimeter integrated navigation method according to claim 4, characterized in that, a k+1 according to t k+1 Time-based velocity estimates v k+1,k The calculation is as follows: ； in, C D The drag coefficient, A For low-orbit satellites, m For low-orbit satellite quality, ρ This refers to the atmospheric density at the location of the low-Earth orbit satellite. lv k+1,k for v k+1,k Size.

6. A terminal device, characterized in that, It includes a processor and a memory; the memory stores computer programs / instructions; the processor executes the computer programs / instructions stored in the memory; the computer programs / instructions are configured to implement the steps of the method according to any one of claims 1 to 5.

Images