Factor graph optimized satellite inertial navigation and ground ranging multiscale fusion positioning method

Abstract

The application discloses a positioning method for satellite inertial navigation and ground ranging multi-scale fusion based on factor graph optimization, comprising the following steps: S1, collecting data in real time by using GNSS, INS and GRT positioning devices; S2, establishing a GNSS and GRT tight combination model based on a factor graph, and solving the tight combination model by using a nonlinear optimization algorithm to obtain a preliminary positioning result of a positioning device; S3, establishing a loose combination model based on a factor graph, fusing the preliminary positioning result obtained in the step S2 with the INS data, and solving the loose combination model by using a nonlinear optimization algorithm to obtain a final positioning result of the positioning device. The application provides a new combined navigation method for accurate navigation and positioning in a satellite insufficient environment, solves a multi-scale fusion difficulty, improves the navigation and positioning precision of the positioning device, and further enhances the positioning robustness in a complex environment.

Description

Technical Field

[0001] This invention belongs to the field of positioning and navigation technology, specifically a multi-scale fusion positioning method based on factor graph optimization of satellite inertial navigation and ground ranging. Background Technology

[0002] With the rapid development of sensor technology and modern communication systems, multi-source autonomous intelligent navigation and positioning has become a core research direction for addressing the challenges of complex environments, especially in scenarios where satellite signals are blocked or interfered with. Achieving high-precision and robust positioning has become the core focus of this field.

[0003] Currently, multi-source fusion navigation and positioning systems are mostly based on satellite navigation, relying on navigation sensors such as inertial navigation and visual navigation as important support. GNSS (Global Navigation Satellite System) can provide accurate all-weather positioning, navigation, and timing services; however, its signals are susceptible to external interference and blockage, leading to decreased positioning accuracy or even complete failure. INS, as an autonomous navigation system, can provide continuously estimated position, velocity, and attitude information, but its navigation errors accumulate over time due to limitations in equipment accuracy and integration principles. Visual and radar-based inertial navigation systems are easily affected by multiple factors such as day / night conditions, clouds, fog, weather, and line-of-sight visibility. Given the performance bottlenecks of single sensors and traditional integrated navigation in complex scenarios, using ground-based ranging (GRT) technologies such as mobile communication, WiFi, and Bluetooth for positioning has gradually become a research hotspot in the navigation field. GRT has abundant signal sources, good coverage, and complete base station facilities; however, these signals are affected by multipath effects and non-line-of-sight propagation during propagation, leading to increased ranging errors and making it difficult to achieve high-precision positioning independently. Therefore, integrating INS, GNSS and GRT signals is of great significance for achieving accurate positioning in complex environments.

[0004] Traditional filter-based loosely coupled GNSS / INS algorithms, based on the assumption of a first-order Markov chain, cannot fully utilize historical observation information and easily degenerate into a pure INS solution mode when satellite signals fail, leading to rapid accumulation and divergence of positioning errors. In contrast, factor graph optimization, as an effective means of solving the maximum a posteriori (MAP) problem, provides a plug-and-play modeling framework. By fully utilizing historical information and relinearization during the iteration process, it achieves more accurate pose estimation than filtering algorithms. Furthermore, factor graphs are naturally compatible with multi-rate, asynchronous, and delayed measurement information from different sensors, while maintaining the sparsity of the MAP problem, significantly improving solution efficiency. Summary of the Invention

[0005] To address the aforementioned problems in existing technologies, the purpose of this invention is to provide a multi-scale fusion positioning method based on factor graph optimization of satellite inertial navigation and ground ranging. A tightly coupled GNSS and GRT model based on factor graphs is established to effectively overcome the multi-scale fusion problem caused by the difference in the magnitude of GNSS and GRT observation data. Simultaneously, to address the multi-rate problem among GNSS, INS, and GRT, a loosely coupled model based on factor graphs is established. This model integrates the advantages of short-term high-precision autonomous navigation of INS with the positioning correction information provided by GNSS and GRT. Through the real-time positioning results provided by GNSS and GRT, the cumulative trend of INS errors is effectively curbed, overcoming the performance limitations of a single sensor in complex environments. This not only improves the navigation and positioning accuracy of the device under positioning but also enhances the stability of positioning, providing strong support for accurate navigation and positioning in environments with insufficient satellite coverage.

[0006] To achieve the above objectives, the technical solution adopted by the present invention is as follows:

[0007] A multi-scale fusion positioning method based on factor graph optimization of satellite inertial navigation and ground ranging includes the following steps: Step S1: Real-time acquisition of GNSS data, INS data, and GRT data of the device to be positioned; Step S2: Establishing a compact combination model of GNSS and GRT based on factor graph, and solving the compact combination model using a nonlinear optimization algorithm to obtain the preliminary positioning result of the device to be positioned; wherein, the compact combination model is used to fuse GNSS pseudorange observations and GRT pseudorange observations; Step S3: Establishing a loose combination model based on factor graph, fusing the preliminary positioning result obtained in step S2 with the INS data, and solving the loose combination model using a nonlinear optimization algorithm to obtain the final positioning result of the device to be positioned.

[0008] As a further improvement to the above technical solution:

[0009] In step S1, the GNSS data includes epoch time, satellite pseudorange measurement value, and the corresponding three-dimensional coordinates of the satellite; the INS data includes epoch time, angular velocity and acceleration of the device to be located; and the GRT data includes epoch time, GRT pseudorange measurement value, and the corresponding three-dimensional coordinates of the base station.

[0010] In step S2, the state variables to be optimized in the GNSS and GRT tightly coupled model based on the factor graph include: the three-dimensional position of the device to be located in the geocentric geofixed coordinate system at epoch k, the GNSS receiver clock deviation, and the GRT receiver clock deviation.

[0011] A GNSS error model is constructed by using GNSS pseudorange measurements and their theoretical predictions, and the covariance matrix of GNSS observations is determined by the satellite carrier-to-noise ratio and satellite elevation angle.

[0012] A GRT error model is constructed by comparing GRT pseudorange measurements with their theoretical predictions, and the covariance matrix of GRT observations is determined by the variance of the GRT ranging error.

[0013] Step S3 specifically includes: transforming the preliminary positioning result of the GNSS and GRT fusion from the geocentric coordinate system to the navigation coordinate system using a rotation matrix to obtain the coordinate-transformed positioning result; performing inertial pre-integration processing on the INS data to obtain the inertial pre-integration formula between adjacent epochs; defining the state variables of the loosely combined model, the state variables including the position, velocity, attitude, gyroscope bias, accelerometer bias, and the coordinate-transformed positioning result of the device to be positioned at epoch k in the navigation coordinate system; constructing a loosely combined error model of the GNSS and GRT positioning result and the INS data based on the factor graph using the inertial pre-integration formula and the coordinate-transformed positioning result; and solving the state variables using a nonlinear optimization algorithm to obtain the final positioning result.

[0014] The nonlinear optimization algorithm is gradient descent, Gauss-Newton method, or Levenberg-Marquardt method, etc.

[0015] In step S1, the ionospheric delay of the GNSS pseudorange measurement is corrected using the Klobuchar model, and the tropospheric delay is compensated using the EGNOS model; the clock deviation of the GRT data is compensated using a time-series prediction model.

[0016] The method is applicable to continuous navigation and positioning of drones, vehicles, or portable devices in complex environments.

[0017] The beneficial effects of this invention are:

[0018] (1) To address the technical challenge of insufficient satellite signals preventing the positioning of devices in complex environments from maintaining high-precision positioning, a multi-scale fusion navigation and positioning method based on factor graph optimization of satellite inertial navigation and ground ranging is proposed. This method effectively overcomes the multi-scale fusion problem caused by the difference in the magnitude of GNSS and GRT data by fusing GNSS and GRT data using factor graphs. Simultaneously, a loosely combined model based on factor graphs is constructed to loosely combine GNSS and GRT positioning results with INS data, resolving the multi-scale and multi-rate problem between GNSS, INS, and GRT. Furthermore, nonlinear optimization algorithms such as gradient descent, Gauss-Newton method, and Levenberg-Marquardt are used to solve the GNSS, INS, and GRT integrated navigation model, resulting in higher-precision positioning results. The real-time positioning results provided by GNSS and GRT effectively curb the accumulation trend of INS errors, overcoming the performance limitations of a single sensor in complex environments.

[0019] (2) By fusing GNSS, INS, and GRT data, continuous high-precision positioning of the device to be positioned can be achieved even when the number of satellites is insufficient. The addition of GRT and INS provides additional sensor information to the device to be positioned, enhancing the navigation and positioning accuracy and stability of the device. Even when the GNSS signal is completely lost, the device to be positioned can achieve continuous and stable positioning through this scheme.

[0020] (3) It effectively solves the technical problem that the device to be positioned cannot maintain high-precision positioning continuously when there are insufficient satellites, improves the accuracy and stability of navigation and positioning, and provides a stable, reliable and accurate navigation and positioning method for the device to be positioned. Attached Figure Description

[0021] Figure 1 This is a schematic diagram of the positioning method of the present invention;

[0022] Figure 2 This is a technical roadmap of the positioning method of the present invention;

[0023] Figure 3 This is a physical image of the device to be located according to an embodiment of the present invention;

[0024] Figure 4 This is a schematic diagram of an experimental scenario according to an embodiment of the present invention;

[0025] Figure 5 This is an experimental result diagram of an embodiment provided by the present invention. Detailed Implementation

[0026] The specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings. It should be understood that the specific embodiments described herein are for illustration and explanation only and are not intended to limit the present invention.

[0027] For ease of description, spatial relative terms such as "above," "on top of," "on the upper surface of," "above," etc., are used herein to describe the spatial positional relationship of a device or feature as shown in the figures to other devices or features. It should be understood that spatial relative terms are intended to encompass different orientations in use or operation beyond the orientation of the device as described in the figures. For example, if the device in the figures were inverted, a device described as "above" or "on top of" other devices or structures would subsequently be positioned as "below" or "under" other devices or structures. Thus, the exemplary term "above" can include both "above" and "below." The device may also be positioned in other different ways (rotated 90 degrees or in other orientations), and the spatial relative descriptions used herein will be interpreted accordingly.

[0028] A multi-scale fusion positioning method based on factor graph optimization of satellite inertial navigation and ground ranging is proposed. This method acquires GNSS and INS data in real time, while simultaneously using ground-based wireless ranging equipment to obtain GRT data. Based on this, a factor graph-based compact combination model of GNSS and GRT is established to overcome the multi-scale fusion problem caused by the difference in the magnitude of GNSS and GRT observation data. Subsequently, to address the multi-rate problem between the positioning results obtained from the compact combination model and INS, a loose combination model based on factor graph is established. Inertial pre-integration is performed on the INS data to obtain the corresponding inertial navigation positioning result at each time step. The preliminary positioning result obtained from the compact combination model of GNSS and GRT is loosely combined with the inertial pre-integration positioning result to obtain continuous and reliable positioning results, effectively improving the navigation and positioning accuracy and stability in complex environments.

[0029] The principle diagram of the positioning method is as follows: Figure 1 As shown, the technology roadmap is as follows: Figure 2 As shown, the specific implementation steps of this method are as follows:

[0030] Step S1: Collect satellite navigation system data (hereinafter referred to as GNSS data) and inertial navigation system data (hereinafter referred to as INS data) of the device to be positioned in real time, and at the same time acquire ground ranging data (hereinafter referred to as GRT data).

[0031] The GNSS data includes epoch times. Satellite pseudorange measurements and satellite position three-dimensional coordinates The INS data includes epoch times. Angular velocity of the device to be positioned and measuring acceleration The GRT data contains epoch times. GRT pseudorange measurement values Three-dimensional coordinates of base station location .

[0032] The device to be located is equipped with an integrated positioning receiver, which has the ability to simultaneously receive and process GNSS satellite signals and GRT ground base station signals. In subsequent models, the GNSS receiver clock bias... Clock bias with GRT receiver This refers to the deviation of the local clock of the same receiver relative to different system times.

[0033] In this embodiment, satellite signals and satellite broadcast ephemeris are acquired in real time, and GNSS pseudorange measurements are obtained based on the measurement information received by the user receiver each time. :

[0034]

[0035] in, , Represents the epoch The total number of visible satellites at that time Represents the epoch The actual pseudorange measurement between the satellite and the receiver at that time, i.e., the satellite pseudorange measurement. This represents the actual distance between the satellite and the receiver. At the speed of light, and These represent the satellite receiver clock skew and the satellite clock skew, respectively. and These represent ionospheric delay and tropospheric delay, respectively. The ionospheric delay is accurately corrected using the Klobuchar model, and the tropospheric delay is effectively compensated using the EGNOS model, thereby obtaining high-precision GNSS pseudorange measurements. This indicates satellite pseudorange measurement noise.

[0036] Ground-based wireless ranging equipment is used to receive GRT data, and the position and clock deviation of the GRT receiver are obtained, including the GRT pseudorange measurement value. Represented as:

[0037]

[0038] In the formula, , Represents the epoch The number of base stations at that time This represents the actual pseudorange measurement value between the base station and the receiver, i.e., the GRT pseudorange measurement value. Indicates base station The actual distance between the receiver and the receiver; At the speed of light, and These represent the clock skew of the GRT receiver and the clock skew of the base station, respectively. A CNN-LSTM-based time-series prediction model is used to estimate and compensate for the base station clock skew online in order to eliminate or reduce the base station clock skew error. This indicates that the receiver measures noise, and the errors are independent of each other.

[0039] Real-time acquisition of inertial navigation system (INS) data yields the measured angular velocity of the device to be positioned (such as a drone). and acceleration The measurement process is subject to various errors. This scheme only considers additive noise and slowly changing zero bias error. The inertial navigation measurement equation is expressed as:

[0040]

[0041]

[0042] In the formula, These are gyroscope measurements. The angular velocity is noise-free. The acceleration measurement value, For noise-free acceleration, Let n be the rotation matrix from the navigation coordinate system (n) to the vehicle coordinate system (b). Represents the gravity vector in the navigation coordinate system, where, and These represent the latitude and altitude of the initial position of the device to be located, respectively. and These represent the zero bias errors of the gyroscope and accelerometer, respectively. and These represent the white noise from the gyroscope and accelerometer, respectively. The angular velocity and acceleration of the inertial navigation system can be obtained from the measurement equations:

[0043]

[0044]

[0045] in, .

[0046] Step S2: Establish a tight combination model of GNSS and GRT based on factor graphs. Define the state variable X, including the current position of the device to be located. Satellite receiver clock bias and GRT receiver clock bias Construct GNSS error models respectively and GRT error model This enables multi-scale fusion of GNSS and GRT data.

[0047] Define the state variables of the GNSS and GRT tightly coupled model. , ,in, Represents the epoch The state variable at time t, For the device to be located in the historical period Location at any given moment This refers to the satellite receiver clock offset (i.e., the GNSS receiver clock offset). To address the clock skew of the GRT receiver, a factor graph-based compact combination model of GNSS and GRT is constructed. :

[0048]

[0049] in, , Represents the total epoch number. and These represent GNSS and GRT at different epochs. Error model at time step; and These represent the GNSS covariance matrix and the GRT covariance matrix, respectively.

[0050] A GNSS error model is constructed based on GNSS pseudorange measurements and model predictions. GNSS pseudorange is a distance measurement obtained by multiplying the time it takes for a satellite signal to travel from the satellite to the receiver by the speed of light. However, due to various error factors, such as satellite clock errors, ionospheric delay, tropospheric delay, multipath effects, and receiver clock errors, there will be deviations between the actual pseudorange measurement and the theoretical model prediction. The GNSS error model can be expressed as:

[0051]

[0052] in, This represents the actual distance between the satellite and the receiver. For the calendar The satellite's three-dimensional coordinates at that time.

[0053] To make more efficient use of GNSS pseudorange measurements for positioning, it is necessary to calculate the GNSS covariance matrix. The GNSS covariance matrix reflects the statistical characteristics of GNSS pseudorange measurement errors, and its magnitude is related to the satellite's carrier-to-noise ratio (CNR) and elevation angle. A higher CNR indicates better signal quality and relatively smaller measurement errors; a higher elevation angle means less influence from atmospheric factors, resulting in smaller measurement errors. Therefore, the GNSS covariance matrix is ​​calculated based on the satellite CNR and elevation angle. .

[0054] Similarly, a GRT error model is constructed based on GRT pseudorange measurements and model predictions. GRT pseudorange is a distance measurement between the device to be located and a ground reference station, obtained through ground-based wireless ranging equipment. It is susceptible to various error factors, such as GRT receiver clock bias and multipath effects caused by the ground environment. Therefore, the GRT error model is defined as:

[0055]

[0056] in, Indicates base station The actual distance between the receiver and the receiver. For the calendar The three-dimensional coordinates of the base station at that time;

[0057] The variance of the GRT ranging error reflects the dispersion of GRT measurements, and its magnitude is related to factors such as the performance of the GRT system and the ground environment. The GRT covariance matrix can be obtained by squaring the variance of the GRT ranging error. .

[0058] The fused position is output using a tight combination model of GNSS and GRT. GNSS receiver clock bias and GRT receiver clock bias .

[0059] Step S3: Establish a loose combination model based on factor graphs. Combine GNSS and GRT positioning results. From the geocentric coordinate system to the navigation coordinate system, define the state variables of the loosely combined model of GNSS and GRT positioning results and INS data based on factor graphs. including location ,speed ,attitude accelerometer zero bias and gyroscope zero bias Pre-integration processing is performed on the INS data to obtain the inertial navigation positioning result at the current integration time, and a loose combination model based on the factor graph is established. and the corresponding covariance matrix The loose combination model was solved using nonlinear optimization algorithms such as gradient descent, Gauss-Newton method, and Levenberg-Marquardt method to obtain the integrated navigation and positioning results.

[0060] The GNSS and GRT positioning results are loosely combined with INS data; specifically, the preliminary positioning results are obtained by tightly combining GNSS and GRT data. By rotation matrix The positioning result is obtained by transforming from the Earth-centered Earth-fixed coordinate system (e-frame) to the navigation coordinate system (n-frame). The rotation matrix is ​​defined as:

[0061]

[0062] In the formula, L represents the longitude of the current positioning result, and B represents the latitude of the current positioning result.

[0063] Pre-integration processing is performed on the INS data, and the derivation is based on the inertial navigation angular velocity and acceleration. Time's up Inertial navigation integral model at any time:

[0064]

[0065] In the formula, The carriers are respectively in Position, velocity, and attitude in the navigation coordinate system at all times. The carriers are respectively in Position, velocity, and attitude in the navigation coordinate system at all times. for The attitude rotation matrix from the inertial navigation coordinate system to the navigation coordinate system at any time; , They represent Time and time, and They are respectively The outputs of the accelerometer and gyroscope are constantly being monitored; For Earth's gravity in the navigation coordinate system, These are the parameters of the noise model for the inertial measurement unit.

[0066] According to the inertial navigation integral model, then from Time's up The inertial pre-integral formula at time t is:

[0067]

[0068] , , They represent from Time's up Pre-integration of position, velocity, and attitude at specific moments. Indicates from Time's up The attitude rotation matrix at time t.

[0069] Define the state variables of a loosely combined model of GNSS and GRT positioning results and INS data based on factor graphs. , Represents the total epoch number. For the calendar The state at time, including the position of the device to be located in the navigation coordinate system at epoch k. ,speed ,attitude gyroscope zero bias accelerometer zero bias and the positioning result after coordinate transformation .

[0070] Construct the following loose combination model based on factor graphs. :

[0071]

[0072] in, Defined as:

[0073]

[0074] In the formula, Indicates the device to be located in the epoch. The state variables at time t, including Instagram Position in the system ,speed ,attitude and gyroscope zero bias accelerometer zero bias Location calculated by GNSS and GRT , Represents the inertial navigation covariance matrix; , , , , Represents the total epoch number. and It is a constant; They represent the position pre-integration, respectively. Velocity pre-integral relative to accelerometer zero bias Jacobian matrix, They represent the position pre-integration, respectively. Speed ​​pre-integral attitude pre-integration relative to gyroscope zero bias Jacobian matrix; yes +1 moment on Instagram Its position in the system, yes +1 moment on Instagram The velocity in the system, yes +1 moment INS from Tie Rotation matrix of the system, Represents the epoch The gyroscope is at zero bias at all times. Represents the epoch The accelerometer is always at zero bias.

[0075] The GNSS and GRT tightly coupled model and the loosely coupled model based on factor graphs employ nonlinear optimization algorithms such as gradient descent, Gauss-Newton method, and Levenberg-Marquardt method for navigation and positioning solutions. Initial state estimates and iteration termination conditions are set. In each iteration, the gradient of the objective function and the Hessian matrix are calculated, and then the state estimates are updated. This process continues until the iteration termination condition is met—that is, the change in the objective function is less than a certain threshold or the number of iterations reaches its maximum value—ultimately yielding continuous positioning results for the integrated navigation.

[0076] The following is an experimental verification of the combined navigation and positioning method of multi-scale fusion of satellite inertial navigation and ground ranging proposed in this invention, which is optimized by factor graph, through a UAV flight experiment.

[0077] like Figure 3 The drone shown is an experimental platform, and the data acquisition equipment and parameters used are shown in Table 1.

[0078] Table 1: GNSS, INS and GRT Data Acquisition Equipment and Parameters

[0079]

[0080] According to the method and model in the invention, an inertial coordinate system, a navigation coordinate system, and a geocentric-ground-fixed coordinate system are defined. The inertial navigation satellite acquisition equipment and the ground wireless ranging equipment are started. During the flight of the UAV, the inertial navigation system outputs acceleration and angular velocity information in the inertial coordinate system at a high frequency of 125Hz. The GNSS receiver outputs pseudorange information and satellite position in the geocentric-ground-fixed coordinate system at a frequency of 1Hz. The ground wireless ranging equipment outputs GRT pseudorange information and the corresponding base station position in the geocentric-ground-fixed coordinate system at a frequency of 1Hz.

[0081] The corresponding GNSS and GRT data are obtained from the GNSS, INS and GRT databases. The preliminary positioning results of the GNSS and GRT compact combination model based on factor graph are solved at the current time. Then the positioning results are transformed from the geocentric geofixed coordinate system to the navigation coordinate system.

[0082] The system obtains relevant inertial navigation data from GNSS, INS, and GRT databases, calculates the positioning result at the current moment through the inertial navigation positioning module, and fuses the inertial navigation positioning result with the GNSS and GRT positioning results through a slack combination model based on factor graphs to finally obtain the accurate positioning result of the UAV.

[0083] To verify the effectiveness of the method proposed in this invention, a UAV flight test was conducted. The test environment was as follows: Figure 4As shown, the second track and field stadium of a university was selected as the flight scenario. This scenario has open space and relatively stable environmental conditions, which is conducive to obtaining accurate positioning data. The UAV first flies in a straight line to the designated position, and then flies around a rectangle to collect data. The rectangular flight trajectory can simulate the complex motion states of the UAV in actual flight, such as turning, acceleration, and deceleration, and comprehensively evaluate the performance of the positioning system under different motion conditions.

[0084] The drone positioning results are as follows Figure 5 As shown, red represents the true positioning value, and blue represents the positioning result based on the fusion of GNSS, INS, and GRT using extended Kalman filter (EKF). Figure 5 The image shows the positioning results (referred to as the comparison method in Chinese), with green representing the fusion positioning results based on the method proposed in this invention. As can be seen from the positioning results image, the positioning results of the method proposed in this invention are closer to the true positioning value, indicating a significant advantage in positioning accuracy. This is because the factor graph-based compact and loose combination models used in this invention can more effectively handle multi-source heterogeneous data, fully considering the correlation and complementarity between different sensor data, and improving the accuracy and reliability of the positioning solution by optimizing the error function. In contrast, the EKF-based fusion method has certain limitations in handling nonlinear problems and multi-source data fusion, resulting in relatively lower positioning accuracy.

[0085] Finally, it is necessary to state that the above embodiments are only used to further illustrate the technical solution of the present invention in detail, and should not be construed as limiting the scope of protection of the present invention. Any non-essential improvements and adjustments made by those skilled in the art based on the above content of the present invention shall fall within the scope of protection of the present invention.

Claims

1. A multi-scale fusion positioning method based on satellite inertial navigation and ground ranging optimized by factor graphs, characterized in that, Includes the following steps: Step S1: Real-time acquisition of GNSS data, INS data, and GRT data of the device to be located; Step S2: Establish a GNSS and GRT compact combination model based on factor graphs, and use a nonlinear optimization algorithm to solve the compact combination model to obtain the preliminary positioning results of the device to be located; wherein, the compact combination model is used to fuse GNSS pseudorange observations and GRT pseudorange observations; Step S3: Establish a loose combination model based on factor graph, fuse the preliminary positioning results obtained in step S2 with the INS data, and use a nonlinear optimization algorithm to solve the loose combination model to obtain the final positioning result of the device to be positioned.

2. The positioning method according to claim 1, characterized in that: In step S1, the GNSS data includes epoch time, satellite pseudorange measurement value, and the corresponding three-dimensional coordinates of the satellite; the INS data includes epoch time, angular velocity and acceleration of the device to be located; and the GRT data includes epoch time, GRT pseudorange measurement value, and the corresponding three-dimensional coordinates of the base station.

3. The positioning method according to claim 1, characterized in that: In step S2, the state variables to be optimized in the GNSS and GRT tight combination model based on the factor graph include: the three-dimensional position of the device to be located in the geocentric geofixed coordinate system at epoch k, the GNSS receiver clock deviation, and the GRT receiver clock deviation. A GNSS error model is constructed by using GNSS pseudorange measurements and their theoretical predictions, and the covariance matrix of GNSS observations is determined by the satellite carrier-to-noise ratio and satellite elevation angle. A GRT error model is constructed by comparing GRT pseudorange measurements with their theoretical predictions, and the covariance matrix of GRT observations is determined by the variance of the GRT ranging error.

4. The positioning method according to claim 1, characterized in that: Step S3 specifically includes: By using a rotation matrix, the preliminary positioning result of the GNSS and GRT fusion is transformed from the geocentric coordinate system to the navigation coordinate system to obtain the positioning result after coordinate transformation. Inertial pre-integration processing is performed on the INS data to obtain the inertial pre-integration formula between adjacent epochs; Define the state variables of the loose combination model, which include the position, velocity, attitude, gyroscope bias, accelerometer bias, and positioning result after coordinate transformation of the device to be positioned in the navigation coordinate system at epoch k. A loose combination error model of GNSS and GRT positioning results and INS data based on factor graphs is constructed using the inertial pre-integration formula and the positioning results after coordinate transformation. The state variables are solved using a nonlinear optimization algorithm to obtain the final positioning result.

5. The positioning method according to claim 1, characterized in that: The nonlinear optimization algorithm is gradient descent, Gauss-Newton method, or Levenberg-Marquardt method.

6. The positioning method according to claim 1, characterized in that: In step S1, the ionospheric delay of the GNSS pseudorange measurement is corrected using the Klobuchar model, and the tropospheric delay is compensated using the EGNOS model; the clock deviation of the GRT data is compensated using a time-series prediction model.

7. The positioning method according to claim 1, characterized in that: The method is applicable to continuous navigation and positioning of drones, vehicles, or portable devices in complex environments.

Images