Positioning method and system based on optical goniometry and astronomic-inertial combined navigation

Abstract

The present application relates to the technical field of integrated navigation, in particular to a positioning method and system based on optical angle measurement and astronomical-inertial integrated navigation. The method first acquires star and navigation satellite images by shooting a star map, and then calculates the attitude of the camera relative to the geocentric inertial coordinate system by identifying and processing the stars, and further determines the coordinates of the navigation satellite in the coordinate system. Through coordinate conversion, the angle information of the navigation satellite relative to the observer in the earth coordinate system is obtained. The optical angle measurement information is fused with the attitude, speed and position information output by the inertial navigation system, the Kalman filter state equation and observation equation are constructed, and finally the real-time high-precision solution and update of the observer's position are realized. The present application combines the accuracy of optical angle measurement, the limited distance reference of astronomical navigation and the autonomous continuity of inertial navigation, effectively improving the navigation and positioning accuracy and reliability under complex environment and satellite denial conditions.

Description

Technical Field

[0001] This invention relates to the field of integrated navigation technology, and in particular to a positioning method and system based on optical angle measurement and astronomical-inertial integrated navigation. Background Technology

[0002] Inertial navigation systems (INS) are characterized by complete autonomy and strong anti-interference capabilities, making them an effective navigation and positioning method even when satellite navigation systems are subject to strong electromagnetic interference (BeiDou / GPS positioning failure). However, their errors diverge over time, failing to meet the requirements for long-endurance navigation and positioning. Celestial navigation, which uses observations of celestial bodies and local horizontal information for positioning, can be used to correct the accumulated position errors of INS. However, traditional celestial navigation systems, due to the infinite distance of stars from Earth, require the use of horizontal information, resulting in lower positioning accuracy.

[0003] In existing technologies, the positioning error of a single inertial navigation system accumulates over time, and astronomical navigation receives parallel starlight, requiring the combination of local horizon information. Furthermore, existing technologies rely on low-Earth orbit (LEO) satellites for optical observation. Compared to navigation satellites, LEO satellites have lower orbital accuracy and poorer predictability, and are more susceptible to Earth's shadow, limiting their observation windows to the early morning and late afternoon, thus failing to meet the requirements for long-endurance, high-precision navigation and positioning.

[0004] Therefore, the present invention provides a positioning method and system based on optical angle measurement and astronomical-inertial integrated navigation to solve the technical problems existing in the prior art. Summary of the Invention

[0005] The purpose of this invention is to provide a positioning method and system based on optical angle measurement and astronomical-inertial integrated navigation. The specific technical solution is as follows:

[0006] The positioning method based on optical angle measurement and astronomical-inertial integrated navigation includes the following steps:

[0007] Step S1: Track navigation satellites based on satellite ephemeris and take star images so that stars and satellites are imaged together in the same star image;

[0008] Step S2: Extract star points from the star map captured in Step S1, match and identify them with the navigation star database, and determine the star coordinates in the geocentric inertial coordinate system. If the number of identified stars exceeds the set threshold, proceed to step S3; otherwise, select a new star map.

[0009] Step S3: Based on the coordinates of the identified stars in the geocentric inertial coordinate system and the star map coordinates... An attitude calculation model was established, and the attitude transformation matrix of the geocentric inertial coordinate system relative to the camera coordinate system was calculated. This allows us to derive the coordinates of the navigation satellite in the geocentric inertial coordinate system.

[0010] Step S4: Based on the coordinates of the navigation satellite in the geocentric inertial coordinate system obtained in step S3, the coordinates of the navigation satellite in the geocentric inertial coordinate system are transformed to the Earth coordinate system through epoch transformation, and the observation angle information of the navigation satellite relative to the observer is calculated.

[0011] Step S5: Based on the observation angle information obtained in step S4, and combined with the attitude, velocity, and position information output by the inertial navigation system, construct the Kalman filter state equation and observation equation.

[0012] Step S6: Based on the Kalman filter state equation and observation equation, the Kalman filter algorithm is used to estimate and update the observer's position information in real time.

[0013] Furthermore, in step S2, the navigation star database includes any one of the Gaia, Hipparcos, and Tycho star catalogs;

[0014] In step S2, star map recognition is performed using the grid method, specifically including:

[0015] ① First determine the main star and mode radius The star pattern is determined by the radius. The system consists of two companion stars, with the remaining companion stars removed as redundant stars.

[0016] ②Shift the star map and reposition the primary star. This places the main star in the center of the field of view;

[0017] ③ Determine the radius of the nearest star And within the radius of the nearest star Outer and mode radius Within, select the distance from the primary star. The nearest companion star is taken as the nearest neighbor star; the field of view is rotated about the line connecting the primary star and the nearest neighbor star, and finally the field of view is divided into... A grid is a feature pattern used to construct stars;

[0018] ④ Match the characteristic patterns of the observed stars with the patterns in the navigation star database to complete star map recognition.

[0019] Furthermore, step S3 specifically involves:

[0020] Satellites within the camera's field of view are imaged onto the detector through the camera lens, with the imaging model approximating a pinhole imaging model; the imaging coordinates of the star points on the detector are obtained through star point extraction. If the coordinates of the camera's optical axis on the detector are focal length is Then the star vector in the camera coordinate system Represented as:

[0021] ;

[0022] Star vectors in the geocentric inertial coordinate system are obtained through star map identification. Ideally, the following relationship exists:

[0023] ;

[0024] Use subscripts Indicates the first star identified within a single frame of the star map One star, ;make and ,but

[0025] ;

[0026] When the number of stars is identified When the value is greater than 2, the camera attitude is solved using the QUEST algorithm, and the attitude transformation matrix between the geocentric inertial coordinate system and the camera coordinate system is obtained. ;

[0027] Vector of navigation satellite in camera coordinate system Substitution The vector of the navigation satellite in the geocentric inertial coordinate system is obtained. This allows us to obtain the coordinates of the navigation satellite in the geocentric inertial coordinate system.

[0028] .

[0029] Furthermore, in step S4, the observers include: a ground carrier with altitude information, a ship moving at sea level, and an aircraft equipped with an altimeter.

[0030] The specific process of epoch conversion includes: from the coordinates of the International Celestial Reference System based on the J2000.0 epoch coordinate system, after correcting for the Earth's proper motion, parallax, offset-precession-nutation and solar gravitational lensing effect, the coordinates of the navigation satellite in the intermediate celestial reference system are obtained; after correcting for the geodynamic parameters and atmospheric refraction, the coordinates of the navigation satellite in the Earth coordinate system are obtained.

[0031] The vector of the navigation satellite in the geocentric inertial coordinate system To transform to Earth coordinates, multiply by the transformation matrix on the left. The vector of the navigation satellite in the Earth coordinate system ;

[0032] Based on the geometric relationship between the observer and the navigation satellite, the azimuth angle of the navigation satellite relative to the observer is calculated. and elevation angle For the first One satellite, azimuth angle and elevation angle Represented as:

[0033] ;

[0034] in: For the first The coordinates of the satellite in the Earth coordinate system These are the observer's coordinates in the Earth coordinate system.

[0035] Furthermore, the bias-precession-nutation matrix Represented as:

[0036] ;

[0037] In the formula, Represents the nutation matrix. Represents the precession matrix, Represents the bias matrix;

[0038] ;

[0039] In the formula, and The celestial polar offset at epoch J2000. The right ascension offset of the J2000 level equatorial coordinate system relative to the geocentric celestial coordinate system; the classical precession is expressed as... , , and ; The obliquity of the ecliptic at epoch J2000; The nutation momentum of the lunar-sun nutation. The nutation momentum of a planet's nutation; , , This represents the rotation matrix about the coordinate axes, with subscripts 1, 2, and 3 corresponding to the values ​​in the Cartesian coordinate system. x, y, z The axis of rotation follows the right-hand rule.

[0040] The geodynamic parameters include polar motion and diurnal aberration;

[0041] Sunday's light difference correction is expressed as:

[0042] ;

[0043] In the formula, For the amplitude of optical aberration, This represents the relative velocity between the observation point and the observed target. This represents the angle between the relative velocity vector and the line of sight. Represents the speed of light;

[0044] The angle between the relative velocity vector and the line of sight is expressed as:

[0045] ;

[0046] In the formula, Represents the velocity in the geocentric inertial coordinate system. Represents the line-of-sight vector in the geocentric inertial coordinate system;

[0047] The unit normal to the plane formed by the relative velocity vector and the view direction Represented as:

[0048] ;

[0049] By using the amplitude and direction of the axial aberration, the transfer matrix of the view direction relative to the true direction can be obtained. :

[0050] ;

[0051] According to the transfer matrix Compensation is applied to the viewing direction:

[0052] ;

[0053] in: This represents the line-of-sight vector in the Earth coordinate system. This represents the true direction obtained after compensating for axial aberration in the downward direction of the Earth coordinate system.

[0054] Atmospheric refraction correction is expressed as:

[0055] ;

[0056] In the formula, The atmospheric refractive index is in its standard form. Zenith distance;

[0057] In non-standard conditions, atmospheric refraction correction value Represented as:

[0058] ;

[0059] In the formula, This is the local air pressure value. This is the local temperature value.

[0060] Furthermore, in step S5, based on the observation angle information obtained in step S4, linearization processing is performed to obtain the error state observation model:

[0061] ;

[0062] in: The observation noise follows a Gaussian distribution and has a non-zero mean; To correct errors in physical quantities; Represents the coefficient matrix. The expression is:

[0063] ;

[0064] in: , ; Indicates the first The distance between the satellite and the observer Indicates the first Satellites and observers in Distance between planes;

[0065] Through position transformation matrix coordinates in the Earth coordinate system Transformation to navigation coordinate system, position transformation matrix for:

[0066] ;

[0067] in: Indicates latitude, Indicates longitude. Indicates altitude, For the Earth's radius, The Earth's eccentricity;

[0068] when When a satellite is observed, the observation equation Represented as:

[0069] ;

[0070] If the latitude and longitude height error is modeled as a random constant, then the one-step state transition matrix is ​​expressed as:

[0071] .

[0072] Furthermore, in step S5, the angular velocity information output by the inertial navigation system is denoted as... The original information is recorded as The attitude, velocity, and position information are denoted as follows: , and ;

[0073] Construct the inertial navigation system update equations, including attitude update equations, velocity update equations, and position update equations, where:

[0074] The attitude update equation is:

[0075] ;

[0076] The velocity update equation is:

[0077] ;

[0078] The position update equation is:

[0079] ;

[0080] in: Represents the rotation matrix from the vehicle frame to the navigation frame; The oblique symmetric matrix representing the rotational angular velocity of the carrier system relative to the navigation system; This represents the rotation of the navigation frame relative to the inertial frame. ,in, The rotational angular velocity of the navigation system caused by the Earth's rotation. , The navigation system rotates due to the curvature of the Earth's surface as the system moves across it. ; Indicates speed under navigation system. , , and These represent the speeds in the east, north, and sky directions, respectively. The latitude of the location of the carrier. The longitude of the location of the carrier. The elevation of the location of the carrier; Represents the gravitational acceleration vector; and These are the radius of curvature of the Earth's meridian and the radius of curvature of its trochanter, respectively. , , , , They represent , , , , The differential.

[0081] Furthermore, in step S5, the Kalman filter state equation and observation equation are constructed as follows:

[0082] The inertial navigation error propagation model is as follows:

[0083] ;

[0084] in: Represents the velocity error matrix. , , and These represent the speed errors in the east, north, and sky directions, respectively. Represents the specific force matrix. , , and These represent the specific force output by the accelerometers in the east, north, and sky directions, respectively. This represents the angular velocity error of the navigation frame relative to the inertial frame. This indicates that the accelerometer has zero bias. , and They represent , and Zero bias of the accelerometer in the direction of travel; This indicates that the gyroscope has zero bias. , and They represent , and The gyroscope has zero bias in direction; , and These represent the errors in latitude, longitude, and altitude, respectively; inertial navigation misalignment angle. , for The differential, for:

[0085] ;

[0086] , and These represent pitch, roll, and yaw misalignment angles, respectively.

[0087] The zero bias of the gyroscope and the zero bias of the accelerometer are modeled as constant values ​​plus white noise, as shown below:

[0088] ;

[0089] in: , They are respectively , The differential;

[0090] Based on the inertial navigation error propagation model and the constant plus white noise model, the Kalman filter state equation is obtained as follows:

[0091] ;

[0092] in: Here is the inertial navigation error state transition matrix; For noise level, This represents the random walk noise of the gyroscope. This indicates the random walk noise of the accelerometer. covariance matrix Recorded as ; This refers to the state quantity of inertial navigation error; for The differential;

[0093] set up For 15-dimensional Kalman filter state variables,

[0094] ;

[0095] Based on the set state variables, the observation equations Rewritten as:

[0096] ;

[0097] in: Represents the measurement matrix. , Indicates measurement noise. covariance matrix Recorded as ;

[0098] The specific expression is:

[0099] ;

[0100] in:

[0101] ;

[0102] ;

[0103] ;

[0104] ;

[0105] ;

[0106] ;

[0107] ;

[0108] ;

[0109] ;

[0110] ;

[0111] in: Indicates the Earth's rotation speed. This represents the rotation matrix from the carrier system to the navigation system.

[0112] Furthermore, in step S6, the Kalman filter algorithm is used to estimate and update the observed azimuth angle in real time. and elevation angle Specifically, this includes time updates and measurement updates;

[0113] The time update is specifically based on: The optimal estimate of the state variables at time step 1 is obtained through the state prediction equation. Predicted values ​​of state variables at time points; based on The optimal estimate of the time error covariance matrix is ​​obtained through the covariance prediction equation. Predicted values ​​of the time error covariance matrix;

[0114] The state prediction equation is expressed as:

[0115] ;

[0116] in, To utilize State prediction at any time Status information at any given moment for The optimal estimate obtained after measurement updates at each time step. for Time's up The state transition matrix at time t;

[0117] The covariance prediction equation is expressed as:

[0118] ;

[0119] in: The predicted covariance matrix, for The error covariance matrix is ​​updated at each time step. Let be the covariance matrix of the system noise;

[0120] The measurement update is specifically as follows: based on Calculate the Kalman filter gain based on the measurement matrix at time point and the predicted value obtained from time update; according to The observations at each time point are used to correct the prior estimates obtained from the time update, resulting in... The posterior optimal estimate at time 1 is obtained, and the posterior covariance is updated simultaneously.

[0121] The filter gain matrix is:

[0122] ;

[0123] in: for The Kalman gain matrix at time t. For the measurement matrix, The measurement noise covariance matrix;

[0124] State estimation equation:

[0125] ;

[0126] in, For the posterior optimal estimate, for Observations at a given time;

[0127] The covariance update equation is:

[0128] ;

[0129] in, It is the identity matrix. Let be the posterior covariance matrix.

[0130] A positioning system based on optical angle measurement and astronomical-inertial integrated navigation for implementing the method described above includes:

[0131] An optical imaging unit is used to track navigation satellites based on satellite ephemeris and capture star maps containing stars and navigation satellites;

[0132] A star image processing unit, connected to the optical imaging unit, is used to extract and identify stars from the star image and calculate the camera attitude and navigation satellite coordinates;

[0133] The navigation calculation unit, connected to the star map processing unit, is used to obtain the coordinates of the navigation satellite in the Earth coordinate system, calculate the observation angle information of the navigation satellite relative to the observer, combine it with the inertial navigation information, construct the Kalman filter state equation and observation equation, execute the Kalman filter algorithm, and use the Kalman filter algorithm to estimate and update the observer's position information in real time.

[0134] The application of the technical solution of the present invention has at least the following beneficial effects:

[0135] The positioning method based on optical angle measurement and astronomical-inertial integrated navigation provided by this invention accurately measures the angle information of navigation satellites through optical angle measurement, and combines the limited-distance observation of astronomical navigation with the autonomous continuity of inertial navigation to form a multi-source information fusion integrated navigation system. This effectively solves the positioning accuracy and continuity problems of single navigation methods in complex environments, and is suitable for long-endurance high-precision navigation and positioning needs in satellite navigation denied environments.

[0136] This invention uses navigation satellites as the observation target. Navigation satellites are closer to Earth, and the direction vector within the observation coordinate system changes with the position of the observation point on the Earth's surface, exhibiting a stronger correlation with the relative position between the observer and the observation target. Compared to low-Earth orbit satellites, navigation satellites have higher orbital accuracy, better predictability, are less affected by Earth's shadow, are not limited by dawn or dusk, and can perform optical observations throughout the night. Compared to stars, the distance between navigation satellites and Earth is fixed, requiring no horizontal information; high-precision positioning can be achieved solely through angle measurement, making it an effective positioning method even under strong electromagnetic interference. Attached Figure Description

[0137] The accompanying drawings, which form part of this invention, are used to provide a further understanding of the invention. The illustrative embodiments of the invention and their descriptions are used to explain the invention and do not constitute an undue limitation of the invention. In the drawings:

[0138] Figure 1 This is a flowchart illustrating the positioning method based on optical angle measurement and astronomical-inertial integrated navigation in an embodiment of the present invention.

[0139] Figure 2 This is a schematic diagram of attitude error estimation in an embodiment of the present invention;

[0140] Figure 3 This is a schematic diagram of speed error estimation in an embodiment of the present invention;

[0141] Figure 4 This is a schematic diagram of position error estimation in an embodiment of the present invention;

[0142] Figure 5 This is a schematic diagram of gyroscope zero bias estimation in an embodiment of the present invention;

[0143] Figure 6 This is a schematic diagram of accelerometer zero bias estimation in an embodiment of the present invention. Detailed Implementation

[0144] The technical solutions in the embodiments of the present invention will be clearly and completely described below. Obviously, the described embodiments are only a part of the embodiments of the present invention, and not all of them. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0145] In the description of this invention, it should be noted that the terms "upper", "lower", "left", "right", "vertical", "horizontal", "inner", "outer", "front", "back", "lateral", "longitudinal", etc., indicate the orientation or positional relationship based on the orientation or positional relationship shown in the accompanying drawings. They are only for the convenience of describing this invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation. Therefore, they should not be construed as limitations on this invention.

[0146] Furthermore, the terms "first," "second," etc., are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of technical features indicated. Thus, a feature defined with "first," "second," etc., may explicitly or implicitly include one or more of that feature. In the description of this invention, unless otherwise stated, "a plurality of" means two or more.

[0147] Example:

[0148] This invention proposes a positioning method based on optical angle measurement and astronomical-inertial integrated navigation. (See [link to relevant documentation]). Figure 1 This includes the following steps:

[0149] Step S1: Track navigation satellites and capture star images based on satellite ephemeris, so that stars and satellites are imaged together in the same star image; satellite ephemeris is obtained by parsing two lines of orbital element files, and star images are captured by a satellite tracking system equipped with a camera.

[0150] Step S2: Extract star points from the star map captured in Step S1, match and identify them with the navigation star database, and identify the stars in the geocentric inertial coordinate system (GCS). Star coordinates in (system) If the number of identified stars exceeds the set threshold, proceed to step S3; otherwise, select a new star map.

[0151] In step S2, the navigation star database includes any one of the Gaia, Hipparcos, and Tycho star catalogs (not limited to these three catalogs);

[0152] Star map identification is performed using a grid method, specifically including:

[0153] ① First determine the main star and mode radius The star pattern is determined by the radius. The system consists of two companion stars, with the remaining companion stars removed as redundant stars.

[0154] ②Shift the star map and reposition the primary star. This places the main star in the center of the field of view;

[0155] ③ Determine the radius of the nearest star And within the radius of the nearest star Outer and mode radius Within, select the distance from the primary star. The nearest companion star is designated as the nearest neighbor; avoid having the nearest neighbor star and the primary star fall on the same grid cell, radius The companion star within will not appear as part of the tectonic star model. The field of view is rotated about the line connecting the primary star and the nearest neighbor star, and finally divided into... A grid is a feature pattern used to construct stars;

[0156] ④ Match the characteristic patterns of the observed stars with the patterns in the navigation star database to complete star map recognition. Specifically:

[0157] Set each grid cell containing stars to 1 and cells without stars to 0, so that the characteristic pattern of the primary star is... grid To indicate, to make Then the first Unit Represented as:

[0158] ;

[0159] Given the observed star pattern and a set of patterns in the navigation star list Star chart recognition is essentially about finding the largest match between observed stars and navigation star catalogs.

[0160] .

[0161] Step S3: Based on the coordinates of the identified stars in the geocentric inertial coordinate system and the star map coordinates... An attitude calculation model was established, and the relative coordinates of the geocentric inertial coordinate system to the camera coordinate system were calculated. c attitude transformation matrix of the system This allows us to derive the coordinates of the navigation satellite in the geocentric inertial coordinate system; specifically:

[0162] Satellites within the camera's field of view are imaged onto the detector through the camera lens, with the imaging model approximating a pinhole camera model; the imaging coordinates of the star points on the detector are obtained through star point extraction. If the coordinates of the camera's optical axis on the detector are focal length is Then the star vector in the camera coordinate system Represented as:

[0163] ;

[0164] Star vectors in the geocentric inertial coordinate system are obtained through star map identification. Ideally, the following relationship exists:

[0165] ;

[0166] Use subscripts Indicates the first star identified within a single frame of the star map One star, ;make and ,but

[0167] ;

[0168] When the number of stars identified When the value is greater than 2, the camera attitude is solved using the QUEST algorithm, and the attitude transformation matrix between the geocentric inertial coordinate system and the camera coordinate system is obtained. Compared with traditional astronauts, this method can directly calculate the local attitude information and simultaneously observe the coordinates of stars in the field of view, thus achieving high-precision attitude calculation.

[0169] Vector of navigation satellite in camera coordinate system Substitution The vector of the navigation satellite in the geocentric inertial coordinate system is obtained. This allows us to obtain the coordinates of the navigation satellite in the geocentric inertial coordinate system.

[0170] .

[0171] Step S4: Based on the coordinates of the navigation satellite in the geocentric inertial coordinate system obtained in step S3, transform the coordinates of the navigation satellite in the geocentric inertial coordinate system to the Earth coordinate system through epoch transformation. Under the system, the observation angle information of the navigation satellite relative to the observer is calculated;

[0172] In step S4, the observers include: a ground carrier with altitude information, a ship moving at sea level, and an aircraft equipped with an altimeter;

[0173] The specific process of epoch conversion includes: from the coordinates of the International Celestial Reference System based on the J2000.0 epoch coordinate system, after correcting for the Earth's proper motion, parallax, offset-precession-nutation, and solar gravitational lensing effect, the coordinates in the intermediate celestial reference system are obtained. Then, after correcting for geodynamic parameters (including polar motion and diurnal aberration) and atmospheric refraction, the coordinates of the navigation satellite in the Earth coordinate system are obtained. In this embodiment, the vector of the navigation satellite in the geocentric inertial coordinate system is... To transform to Earth coordinates, multiply by the transformation matrix on the left. The vector of the navigation satellite in the Earth coordinate system ;

[0174] In this embodiment, the bias-precession-nutation matrix Represented as:

[0175] ;

[0176] In the formula, Represents the nutation matrix. Represents the precession matrix, Represents the bias matrix;

[0177] ;

[0178] In the formula, and The celestial polar offset at epoch J2000. The right ascension offset of the J2000 level equatorial coordinate system relative to the geocentric celestial coordinate system; the classical precession is expressed as... , , and ; The obliquity of the ecliptic at epoch J2000; The nutation momentum of the lunar-sun nutation. The nutation momentum of a planet's nutation; , , This represents the rotation matrix about the coordinate axes, with subscripts 1, 2, and 3 corresponding to the values ​​in the Cartesian coordinate system. x, y, z The axis of rotation follows the right-hand rule.

[0179] The geodynamic parameters include polar motion and diurnal aberration;

[0180] Sunday's light difference correction is expressed as:

[0181] ;

[0182] In the formula, For the amplitude of optical aberration, This represents the relative velocity between the observation point and the observed target. This represents the angle between the relative velocity vector and the line of sight. Represents the speed of light;

[0183] The angle between the relative velocity vector and the line of sight is expressed as:

[0184] ;

[0185] In the formula, Represents the velocity in the geocentric inertial coordinate system. Represents the line-of-sight vector in the geocentric inertial coordinate system;

[0186] The unit normal to the plane formed by the relative velocity vector and the view direction Represented as:

[0187] ;

[0188] By using the amplitude and direction of the axial aberration, the transfer matrix of the view direction relative to the true direction can be obtained. :

[0189] ;

[0190] According to the transfer matrix Compensation is applied to the viewing direction:

[0191] ;

[0192] in: This represents the line-of-sight vector in the Earth coordinate system. This represents the true direction obtained after compensating for axial aberration in the downward direction of the Earth coordinate system.

[0193] Atmospheric refraction correction is expressed as:

[0194] ;

[0195] In the formula, The atmospheric refractive index is in its standard form. Zenith distance;

[0196] In non-standard conditions, atmospheric refraction correction value Represented as:

[0197] ;

[0198] In the formula, This is the local air pressure value. This is the local temperature value.

[0199] Based on the geometric relationship between the observer and the navigation satellite, the azimuth angle of the navigation satellite relative to the observer is calculated. and elevation angle For the first One satellite, azimuth angle and elevation angle Represented as:

[0200] ;

[0201] in: For the first The coordinates of the satellite in the Earth coordinate system These are the observer's coordinates in the Earth coordinate system.

[0202] Step S5: Based on the observation angle information obtained in step S4, and combined with the attitude, velocity, and position information output by the inertial navigation system, construct the Kalman filter state equation and observation equation; specifically:

[0203] Based on the observation angle information obtained in step S4, a linearization process is performed to obtain the error state observation model:

[0204] ;

[0205] in: The observation noise follows a Gaussian distribution and has a non-zero mean; To correct errors in physical quantities; Represents the coefficient matrix. The expression is:

[0206] ;

[0207] in: , ; Indicates the first The distance between the satellite and the observer Indicates the first Satellites and observers in Distance between planes;

[0208] Through position transformation matrix coordinates in the Earth coordinate system Transform to navigation coordinate system ( (System), position transformation matrix for:

[0209] ;

[0210] in: Indicates latitude, Indicates longitude. Indicates altitude, For the Earth's radius, Earth's eccentricity;

[0211] when When a satellite is observed, the observation equation Represented as:

[0212] ;

[0213] If the latitude and longitude height error is modeled as a random constant, then the one-step state transition matrix is ​​expressed as:

[0214] .

[0215] In step S5, the angular velocity information output by the inertial navigation system is denoted as... The original information is recorded as The attitude, velocity, and position information are denoted as follows: , and ;

[0216] Construct the inertial navigation system update equations, including attitude update equations, velocity update equations, and position update equations, where:

[0217] The attitude update equation is:

[0218] ;

[0219] The velocity update equation is:

[0220] ;

[0221] The position update equation is:

[0222] ;

[0223] in: Represents the rotation matrix from the vehicle frame to the navigation frame; The oblique symmetric matrix representing the rotational angular velocity of the carrier system relative to the navigation system; This indicates the rotation of the navigation frame relative to the inertial frame. ,in, The rotational angular velocity of the navigation system caused by the Earth's rotation. , The navigation system rotates due to the curvature of the Earth's surface as the system moves across it. ; Indicates speed under navigation system. , , and These represent the speeds in the east, north, and sky directions, respectively; location information includes latitude, longitude, and altitude. The latitude of the location of the carrier. The longitude of the location of the carrier. The elevation of the location of the carrier; Represents the gravitational acceleration vector; and These are the radius of curvature of the Earth's meridian and the radius of curvature of its trochanter, respectively. , , , , They represent , , , , The differential.

[0224] In step S5, the Kalman filter state equation and observation equation are constructed as follows:

[0225] The inertial navigation error propagation model is as follows:

[0226] ;

[0227] in: Represents the velocity error matrix. , , and These represent the speed errors in the east, north, and sky directions, respectively. Represents the specific force matrix. , , and These represent the specific force output by the accelerometers in the east, north, and sky directions, respectively. This represents the angular velocity error of the navigation frame relative to the inertial frame. This indicates that the accelerometer has zero bias. , and They represent , and Zero bias of the accelerometer in the direction of travel; This indicates that the gyroscope has zero bias. , and They represent , and The gyroscope has zero bias in direction; , and These represent the errors in latitude, longitude, and altitude, respectively; inertial navigation misalignment angle. , for The differential, for:

[0228] ;

[0229] , and These represent pitch, roll, and heading misalignment angles, respectively.

[0230] The zero bias of the gyroscope and the zero bias of the accelerometer are modeled as constant values ​​plus white noise, as shown below:

[0231] ;

[0232] in: , They are respectively , The differential;

[0233] Based on the inertial navigation error propagation model and the constant plus white noise model, the Kalman filter state equation is obtained as follows:

[0234] ;

[0235] in: Here is the inertial navigation error state transition matrix; For noise level, This represents the random walk noise of the gyroscope. This indicates the random walk noise of the accelerometer. covariance matrix Recorded as ; This refers to the state quantity of inertial navigation error; for The differential;

[0236] set up For 15-dimensional Kalman filter state variables,

[0237] ;

[0238] Based on the set state variables, the observation equations Rewritten as:

[0239] ;

[0240] in: Represents the measurement matrix. , Indicates measurement noise. covariance matrix Recorded as ;

[0241] The specific expression is:

[0242] ;

[0243] in:

[0244] ;

[0245] ;

[0246] ;

[0247] ;

[0248] ;

[0249] ;

[0250] ;

[0251] ;

[0252] ;

[0253] ;

[0254] in: Indicates the Earth's rotation speed. This represents the rotation matrix from the carrier system to the navigation system.

[0255] Step S6: Based on the Kalman filter state equation and observation equation, the Kalman filter algorithm is used to estimate and update the observer's position information in real time; specifically:

[0256] The Kalman filter algorithm is used to estimate and update the observed azimuth in real time. and elevation angle Specifically, this includes time updates and measurement updates;

[0257] The time update is specifically based on: The optimal estimate of the state variables at time step 1 is obtained through the state prediction equation. Predicted values ​​of state variables at time points; based on The optimal estimate of the time error covariance matrix is ​​obtained through the covariance prediction equation. Predicted values ​​of the time error covariance matrix;

[0258] The state prediction equation is expressed as:

[0259] ;

[0260] in, To utilize State prediction at any time Status information at any given moment for The optimal estimate obtained after measurement updates at each time step. for Time's up The state transition matrix at time t;

[0261] The covariance prediction equation is expressed as:

[0262] ;

[0263] in: The predicted covariance matrix, for The error covariance matrix is ​​updated at each time step. Let be the covariance matrix of the system noise;

[0264] The measurement update is specifically as follows: based on Calculate the Kalman filter gain based on the measurement matrix at time point and the predicted value obtained from time update; according to The observations at each time point are used to correct the prior estimates obtained from the time update, resulting in... The posterior optimal estimate at time 1 is obtained, and the posterior covariance is updated simultaneously.

[0265] The filter gain matrix is:

[0266] ;

[0267] in: for The Kalman gain matrix at time t. For the measurement matrix, The measurement noise covariance matrix;

[0268] State estimation equation:

[0269] ;

[0270] in, For the posterior optimal estimate, for Observations at a given time;

[0271] The covariance update equation is:

[0272] ;

[0273] in, It is the identity matrix. Let be the posterior covariance matrix.

[0274] To verify the performance and advantages of the method of the present invention, a simulation experiment was conducted based on data from a certain sea trial. The simulation conditions are set as shown in Table 1 below:

[0275] Table 1 Simulation conditions

[0276]

[0277] The simulation uses the actual ship's sea trial trajectory to generate the observer's real motion trajectory, and superimposes a 100 m initial position error as the simulation trajectory. Using a simulated satellite as the observation satellite, the observation angles (azimuth and elevation) at each position point are generated based on the simulated trajectory. The satellite orbit radial accuracy is set according to the orbit determination accuracy of geostationary orbit satellites, with an orbital altitude of 20,000 km selected. Specific coordinates are shown in Table 2.

[0278] Table 2 Simulated Satellite Positions

[0279]

[0280] Figure 2-6 It demonstrates the accuracy estimation of attitude, velocity, position, gyroscope bias, and accelerometer bias in an astronomical-inertial integrated navigation system.

[0281] Among them, by Figure 2 As can be seen, the final attitude determination error can converge to near 0, demonstrating good attitude determination performance.

[0282] Depend on Figure 3 and Figure 4 The simulation results over 30 hours show that the eastward velocity error is only 0.127 m / s, the northward velocity error is only 0.159 m / s, the eastward positioning error is 17.23 m, and the northward positioning error is 139.25 m, demonstrating good velocity and positioning accuracy, which can meet the positioning requirements.

[0283] Depend on Figure 5 and Figure 6 It can be seen that the gyroscope zero bias eventually converges to 0.003° / h, which is basically consistent with the initial error setting. The accelerometer horizontal zero bias estimate is 35-42 mGal, which shows good zero bias estimation performance and can accurately eliminate zero bias error in practical applications.

[0284] The above experimental results demonstrate the effectiveness of the positioning method based on optical angle measurement and astronomical-inertial integrated navigation proposed in this invention, providing a combined navigation positioning method with good positioning accuracy for all-night observation.

[0285] This invention also provides a positioning system based on optical angle measurement and astronomical-inertial integrated navigation, used to implement the positioning method described above, comprising:

[0286] An optical imaging unit is used to track navigation satellites based on satellite ephemeris and capture star maps containing stars and navigation satellites;

[0287] A star image processing unit, connected to the optical imaging unit, is used to extract and identify stars from the star image and calculate the camera attitude and navigation satellite coordinates;

[0288] The navigation calculation unit, connected to the star map processing unit, is used to obtain the coordinates of the navigation satellite in the Earth coordinate system, calculate the observation angle information of the navigation satellite relative to the observer, combine it with the inertial navigation information, construct the Kalman filter state equation and observation equation, execute the Kalman filter algorithm, and use the Kalman filter algorithm to estimate and update the observer's position information in real time.

[0289] The present invention also provides an electronic device corresponding to the above embodiments. The electronic device may be a processing device for a client, such as a mobile phone, a laptop, a tablet computer, a desktop computer, etc., to execute the methods of the above embodiments.

[0290] The electronic device of this embodiment includes a memory, a processor, and a computer program stored in the memory; the processor executes the computer program in the memory to implement the steps of the method described in the above embodiment.

[0291] In some implementations, the memory may be high-speed random access memory (RAM), and may also include nonvolatile memory, such as at least one disk storage.

[0292] In other implementations, the processor can be any type of general-purpose processor, such as a central processing unit (CPU) or a digital signal processor (DSP), and there is no limitation here.

[0293] The present invention also provides a readable storage medium corresponding to the above embodiments, wherein a computer program / instructions are stored thereon. When the computer program / instructions are executed by a processor, they implement the steps of the methods described in the above embodiments.

[0294] A computer-readable storage medium can be a tangible device that holds and stores instructions for use by an instruction execution device. A computer-readable storage medium can be, for example, but not limited to, an electrical storage device, a magnetic storage device, an optical storage device, an electromagnetic storage device, a semiconductor storage device, or any combination thereof.

[0295] Those skilled in the art will understand that embodiments of this application can be provided as methods, systems, or computer program products. Therefore, this application can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, this application can take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code. The solutions in the embodiments of this application can be implemented in various computer languages, such as the object-oriented programming language Java and the interpreted scripting language JavaScript.

[0296] This invention is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart illustrations and / or block diagrams. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0297] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0298] The above description is only a preferred embodiment of the present invention and does not limit the scope of the present invention. All equivalent structural transformations made using the present invention under the inventive concept of the present invention, or direct / indirect applications in other related technical fields, are included within the protection scope of the present invention.

Claims

1. A positioning method based on optical angle measurement and astronomical-inertial integrated navigation, characterized in that, Includes the following steps: Step S1: Track navigation satellites based on satellite ephemeris and take star images so that stars and satellites are imaged together in the same star image; Step S2: Extract star points from the star map captured in Step S1, match and identify them with the navigation star database, and determine the star coordinates in the geocentric inertial coordinate system. ; If the number of identified stars exceeds the set threshold, proceed to step S3; otherwise, select a new star map. Step S3: Based on the coordinates of the identified stars in the geocentric inertial coordinate system and the star map coordinates... An attitude calculation model was established, and the attitude transformation matrix of the geocentric inertial coordinate system relative to the camera coordinate system was calculated. This allows us to derive the coordinates of the navigation satellite in the geocentric inertial coordinate system. Step S4: Based on the coordinates of the navigation satellite in the geocentric inertial coordinate system obtained in step S3, the coordinates of the navigation satellite in the geocentric inertial coordinate system are transformed to the Earth coordinate system through epoch transformation, and the observation angle information of the navigation satellite relative to the observer is calculated. Step S5: Linearize the observation angle information obtained in step S4 to obtain the error state observation model. Combine the angular velocity information, specific force information, attitude, velocity and position information output by the inertial navigation system to construct the Kalman filter state equation and observation equation. Step S6: Based on the Kalman filter state equation and observation equation, the Kalman filter algorithm is used to estimate and update the observer's position information in real time; Step S3 is as follows: Satellites within the camera's field of view are imaged onto the detector through the camera lens, with the imaging model approximating a pinhole imaging model; the imaging coordinates of the star points on the detector are obtained through star point extraction. If the coordinates of the camera's optical axis on the detector are focal length is Then the star vector in the camera coordinate system Represented as: ； Star vectors in the geocentric inertial coordinate system are obtained through star map identification. Ideally, the following relationship exists: ； Use subscripts Indicates the first star identified within a single frame of the star map One star, ;make and ,but ； When the number of stars is identified When the value is greater than 2, the camera attitude is solved using the QUEST algorithm, and the attitude transformation matrix between the geocentric inertial coordinate system and the camera coordinate system is obtained. ; Vector of navigation satellite in camera coordinate system Substitution The vector of the navigation satellite in the geocentric inertial coordinate system is obtained. This allows us to obtain the coordinates of the navigation satellite in the geocentric inertial coordinate system. 。 2. The positioning method based on optical angle measurement and astronomical-inertial integrated navigation as described in claim 1, characterized in that, In step S2, the navigation star database includes any one of the Gaia star catalog, the Hipparcos star catalog, and the Tycho star catalog; Step S2 specifically includes: ① First determine the main star and mode radius The star pattern is determined by the radius. The system consists of two companion stars, with the remaining companion stars removed as redundant stars. ②Shift the star map and reposition the primary star. This places the main star in the center of the field of view; ③ Determine the radius of the nearest star And within the radius of the nearest star Outer and mode radius Within, select the distance from the primary star. The nearest companion star is taken as the nearest neighbor star; the field of view is rotated about the line connecting the primary star and the nearest neighbor star, and finally the field of view is divided into... A grid is a feature pattern used to construct stars; ④ Match the characteristic patterns of the observed stars with the patterns in the navigation star database to complete star map recognition.

3. The positioning method based on optical angle measurement and astronomical-inertial integrated navigation as described in claim 2, characterized in that, In step S4, the observers include: a ground carrier with altitude information, a ship moving at sea level, and an aircraft equipped with an altimeter; The specific process of epoch conversion includes: from the coordinates of the International Celestial Reference System based on the J2000.0 epoch coordinate system, after correcting for the Earth's proper motion, parallax, offset-precession-nutation and solar gravitational lensing effect, the coordinates of the navigation satellite in the intermediate celestial reference system are obtained; after correcting for the geodynamic parameters and atmospheric refraction, the coordinates of the navigation satellite in the Earth coordinate system are obtained. The vector of the navigation satellite in the geocentric inertial coordinate system To transform to Earth coordinates, multiply by the transformation matrix on the left. The vector of the navigation satellite in the Earth coordinate system ; Based on the geometric relationship between the observer and the navigation satellite, the azimuth angle of the navigation satellite relative to the observer is calculated. and elevation angle For the first One satellite, azimuth angle and elevation angle Represented as: ； in: For the first The coordinates of the satellite in the Earth coordinate system These are the observer's coordinates in the Earth coordinate system.

4. The positioning method based on optical angle measurement and astronomical-inertial integrated navigation as described in claim 3, characterized in that, Bias-precession-nutation matrix Represented as: ； In the formula, Represents the nutation matrix. Represents the precession matrix, Represents the bias matrix; ； In the formula, and The celestial polar offset at epoch J2000. The right ascension offset of the J2000 level equatorial coordinate system relative to the geocentric celestial coordinate system; the classical precession is expressed as... , , and ; The obliquity of the ecliptic at epoch J2000; The nutation momentum of the lunar-sun nutation. The nutation momentum of a planet's nutation; , , This represents the rotation matrix about the coordinate axes, with subscripts 1, 2, and 3 corresponding to the values ​​in the Cartesian coordinate system. x, y, z The axis of rotation follows the right-hand rule. The geodynamic parameters include polar motion and diurnal aberration; Sunday's light difference correction is expressed as: ； In the formula, For the amplitude of optical aberration, This represents the relative velocity between the observation point and the observed target. This represents the angle between the relative velocity vector and the line of sight. Represents the speed of light; The angle between the relative velocity vector and the line of sight is expressed as: ； In the formula, Represents the velocity in the geocentric inertial coordinate system. Represents the line-of-sight vector in the geocentric inertial coordinate system; The unit normal to the plane formed by the relative velocity vector and the view direction Represented as: ； By using the amplitude and direction of the axial aberration, the transfer matrix of the view direction relative to the true direction can be obtained. : ； According to the transfer matrix Compensation is applied to the viewing direction: ； in: This represents the line-of-sight vector in the Earth coordinate system. This represents the true direction obtained after compensating for axial aberration in the downward direction of the Earth coordinate system. Atmospheric refraction correction is expressed as: ； In the formula, The atmospheric refractive index is in its standard form. Zenith distance; In non-standard conditions, atmospheric refraction correction value Represented as: ； In the formula, This is the local air pressure value. This is the local temperature value.

5. The positioning method based on optical angle measurement and astronomical-inertial integrated navigation as described in claim 4, characterized in that, In step S5, based on the observation angle information obtained in step S4, linearization processing is performed to obtain the error state observation model: ； in: The observation noise follows a Gaussian distribution and has a non-zero mean; To correct errors in physical quantities; Represents the coefficient matrix. The expression is: ； in: , ; Indicates the first l The distance between the satellite and the observer Indicates the first l Satellites and observers in xy Distance between planes; Through position transformation matrix coordinates in the Earth coordinate system Transformation to navigation coordinate system, position transformation matrix for: ； in: Indicates latitude, Indicates longitude. Indicates altitude, For the Earth's radius, The Earth's eccentricity; when When a satellite is observed, the observation equation Represented as: ； If the latitude and longitude height error is modeled as a random constant, then the one-step state transition matrix is ​​expressed as: 。 6. The positioning method based on optical angle measurement and astronomical-inertial integrated navigation as described in claim 5, characterized in that, In step S5, the angular velocity information output by the inertial navigation system is denoted as... The original information is recorded as The attitude, velocity, and position information are denoted as follows: , and ; Construct the inertial navigation system update equations, including attitude update equations, velocity update equations, and position update equations, where: The attitude update equation is: ； The velocity update equation is: ； The position update equation is: ； in: Represents the rotation matrix from the vehicle frame to the navigation frame; The oblique symmetric matrix representing the rotational angular velocity of the carrier system relative to the navigation system; This represents the rotation of the navigation frame relative to the inertial frame. ,in, The rotational angular velocity of the navigation system caused by the Earth's rotation. , The navigation system rotates due to the curvature of the Earth's surface as the system moves across it. ; Indicates speed under navigation system. , , and These represent the speeds in the east, north, and sky directions, respectively. The latitude of the location of the carrier. The longitude of the location of the carrier. The elevation of the location of the carrier; Represents the gravitational acceleration vector; and These are the radius of curvature of the Earth's meridian and the radius of curvature of its trochanter, respectively. , , , , They represent , , , , The differential.

7. The positioning method based on optical angle measurement and astronomical-inertial integrated navigation as described in claim 6, characterized in that, In step S5, the Kalman filter state equation and observation equation are constructed as follows: The inertial navigation error propagation model is as follows: ； in: Represents the velocity error matrix. , , and These represent the speed errors in the east, north, and sky directions, respectively. Represents the specific force matrix. , , and These represent the specific force output by the accelerometers in the east, north, and sky directions, respectively. This represents the angular velocity error of the navigation frame relative to the inertial frame. This indicates that the accelerometer has zero bias. , and They represent , and Zero bias of the accelerometer in the direction of travel; This indicates that the gyroscope has zero bias. , and They represent , and The gyroscope has zero bias in direction; , and These represent the errors in latitude, longitude, and altitude, respectively; inertial navigation misalignment angle. , for The differential, for: ； , and These represent pitch, roll, and yaw misalignment angles, respectively. The zero bias of the gyroscope and the zero bias of the accelerometer are modeled as constant values ​​plus white noise, as shown below: ； in: , They are respectively , The differential; Based on the inertial navigation error propagation model and the constant plus white noise model, the Kalman filter state equation is obtained as follows: ； in: Here is the inertial navigation error state transition matrix; For noise level, This represents the random walk noise of the gyroscope. This indicates the random walk noise of the accelerometer. covariance matrix Recorded as ; This refers to the state quantity of inertial navigation error; for The differential; set up For 15-dimensional Kalman filter state variables, ； Based on the set state variables, the observation equations Rewritten as: ； in: Represents the measurement matrix. , Indicates measurement noise. covariance matrix Recorded as ; The specific expression is: ； in: ； ； ； ； ； ； ； ； ； ； in: Indicates the Earth's rotation speed. This represents the rotation matrix from the carrier system to the navigation system.

8. A positioning system for implementing the method described in any one of claims 1-7, based on optical angle measurement and astronomical-inertial integrated navigation, characterized in that, include: An optical imaging unit is used to track navigation satellites based on satellite ephemeris and capture star maps containing stars and navigation satellites; A star image processing unit, connected to the optical imaging unit, is used to extract and identify stars from the star image and calculate the camera attitude and navigation satellite coordinates; The navigation calculation unit, connected to the star map processing unit, is used to obtain the coordinates of the navigation satellite in the Earth coordinate system, calculate the observation angle information of the navigation satellite relative to the observer, combine it with the inertial navigation information, construct the Kalman filter state equation and observation equation, execute the Kalman filter algorithm, and use the Kalman filter algorithm to estimate and update the observer's position information in real time.

Images