A GNSS anti-rolling heading stabilization measurement method based on cross-epoch difference for ships

Abstract

The application discloses a kind of ship-based GNSS anti-tossing heading stability measurement methods based on cross-epoch difference, comprising: constructing double-layer sliding cache queue;Adjacent epochs are executed to carrier phase difference solution short-interval reference heading;Filtering history epoch that meets time or distance threshold as reference epoch, after ephemeris orbit error and atmospheric delay compensation, cross-epoch difference is solved long-interval coarse heading;With short-interval reference heading as the check reference, long-interval coarse heading is executed to heading autonomous integrity detection and outputs stable heading.Through long-short double-layer epoch difference architecture and heading autonomous integrity detection, high stability output and abnormal autonomous identification are realized in single-antenna condition under the condition of tossing working condition heading, the problem that the existing single-antenna direction finding scheme is solved when low-speed tossing heading is violently shaken and lacks integrity monitoring.

Description

Technical Field

[0001] This invention relates to the field of radio direction finding technology, and in particular to a method for measuring the anti-turbulence heading stability of marine GNSS based on cross-epoch difference. Background Technology

[0002] During maritime navigation, heading is one of the core parameters for ensuring navigational safety, trajectory planning, and maneuvering decisions. The accuracy and stability of heading measurement directly affect the reliability of the ship's automated navigation system. Currently, the Global Navigation Satellite System (GNSS) is one of the main means of acquiring ship heading information. Among them, heading measurement schemes based on single-antenna GNSS receivers are widely used in various small and medium-sized civilian and commercial vessels due to their simple hardware structure and low cost. However, under complex operating conditions such as rough seas and low-speed ship navigation, existing single-antenna GNSS heading measurement technology faces prominent problems such as severe heading output jitter and insufficient reliability, which seriously restricts the ship's operational capabilities in high-precision scenarios such as low-speed berthing, fixed-point operations, and channel fine-tuning.

[0003] Currently, mainstream marine GNSS heading measurement methods are mainly divided into two categories. The first category is the Doppler instantaneous velocity vector-based heading calculation scheme. This scheme relies on Doppler observations at the current epoch to calculate the instantaneous eastward and northward velocities of the vessel, and then calculates the real-time heading angle through the angle between the velocity vectors. However, under conditions of low-speed ship navigation and rough seas, the instantaneous velocity disturbance errors caused by irregular swaying and vertical undulations are numerically much larger than the ship's actual speed. The noise in the velocity vector is significantly amplified, leading to continuous and large jumps in the heading angle, with fluctuations reaching ±8° to ±15°. Since Doppler velocity measurement only characterizes the instantaneous motion trend of the vessel at the current epoch, its time-domain smoothing methods are limited, and it cannot suppress heading noise caused by turbulence at the mechanistic level. Therefore, this scheme is basically unusable in high-precision scenarios such as low-speed berthing and fixed-point operations. The second category is the dual-antenna GNSS baseline orientation scheme. This scheme uses a baseline formed by two antennas to perform carrier phase differential calculation, directly calculating the vessel's heading, resulting in relatively high heading stability. However, this solution has obvious shortcomings in application: the hardware cost increases exponentially, requiring two independent radio frequency channels and high requirements for the wiring and installation space of the dual antennas; the heading accuracy is highly dependent on the baseline length, and small ships or simple work vessels often cannot reserve a sufficiently long baseline installation position; in addition, the deformation of the hull and the torsion of the antenna under the action of wind and waves will continuously introduce baseline system errors, and the subsequent calibration and maintenance process is complicated, resulting in high overall implementation costs.

[0004] To improve the jitter problem in single-antenna GNSS heading measurements, existing technologies have proposed improved schemes based on filtering smoothing or multi-sensor fusion. Filtering smoothing algorithms smooth the heading output using methods such as time-domain low-pass filtering or Kalman filtering; however, these methods inevitably introduce a time lag in the heading output, making it difficult to meet the real-time heading response speed requirements of ship maneuvers. Multi-sensor fusion schemes introduce auxiliary sensors such as inertial measurement units (IMUs) or electronic compasses to combine with GNSS for navigation. While this can improve heading stability to some extent, it increases hardware costs and system complexity, and the inherent cumulative drift error of inertial devices cannot fundamentally eliminate the disturbance of sea surface turbulence to GNSS heading measurements.

[0005] In addition, a few existing technologies employ cross-epoch differential methods for heading measurement. These methods utilize the carrier phase difference between the current epoch and a historical epoch to solve for the relative displacement vector, thereby calculating the heading. While these methods leverage long baseline distances to dilute instantaneous disturbance errors and offer some heading stabilization, key technological gaps remain. Firstly, existing methods rely solely on long-interval cross-epoch differentials to output a single heading result, lacking an independent heading verification benchmark and the ability to monitor heading integrity. They cannot autonomously identify heading anomalies caused by cycle slips, multipath errors, or residual orbital errors, making it difficult to guarantee the reliability of the output results. Secondly, when the time interval is long (e.g., exceeding 5 seconds), slowly varying errors such as orbital deviations from satellite broadcast ephemeris, satellite clock errors, and ionospheric and tropospheric delays accumulate continuously with the length of the time interval. Errors that are negligible in short-interval differentials can form systematic relative displacement deviations in long-interval differentials, directly degrading the accuracy of the heading calculation. Existing cross-epoch difference schemes lack a complete ephemeris orbit and atmospheric error compensation process for long time difference scenarios, which limits the applicable time range and accuracy limit of cross-epoch difference. Summary of the Invention

[0006] The purpose of this invention is to provide a method for measuring the heading stability of marine GNSS under turbulence based on cross-epoch differential. By constructing a long and short double-layer epoch differential architecture and using the short-interval reference heading as the verification benchmark, the method performs autonomous heading integrity detection on the long-interval coarse heading. Under the condition of low cost with a single antenna, it achieves high stability output of heading and autonomous identification of anomalies under turbulence conditions. This solves the problem that existing single-antenna GNSS direction finding schemes suffer from severe heading jitter and lack the ability to monitor the integrity of the heading layer when the ship is under low-speed turbulence.

[0007] To address the aforementioned technical problems, a first aspect of this invention provides a method for measuring the anti-turbulence heading stability of marine GNSS based on cross-epoch difference, comprising the following steps: S1. Real-time acquisition of raw observation data from shipborne GNSS receivers to construct a two-layer sliding buffer queue. The two-layer sliding buffer queue includes a long-interval buffer area for caching observation data of multiple historical epochs, and a short-interval buffer area for caching observation data of the current epoch and the previous adjacent epoch. S2, perform carrier phase difference between the current epoch and the previous adjacent epoch in the short interval buffer to calculate the short interval reference heading; S3. Select historical epochs that meet at least one of the preset time threshold and preset horizontal distance threshold from the long interval buffer as reference epochs. After performing ephemeris orbital error compensation and atmospheric delay compensation on the current epoch and the reference epoch, perform cross-epoch carrier phase differential to calculate the long interval coarse heading. S4, using the short-interval reference heading as a verification benchmark, perform heading autonomous integrity detection on the long-interval coarse heading to obtain a stable heading output for the current epoch; the heading autonomous integrity detection includes: calculating the heading difference between the short-interval reference heading and the long-interval coarse heading; if the heading difference exceeds a preset integrity threshold, then determine that the long-interval coarse heading is abnormal and output the short-interval reference heading; otherwise, output the long-interval coarse heading.

[0008] Further, the calculation of the short-interval reference heading in step S2 includes: S21, Perform epoch-to-epoch difference on the carrier phase observations of the current epoch and the previous adjacent epoch to eliminate carrier integer ambiguity and solve for the relative position increment in the short-interval geocentric geofixed coordinate system; S22, the relative position increment in the short-interval geocentric-ground-fixed coordinate system is transformed to the ship's local northeast-sky coordinate system, the vertical undulation displacement is separated, the horizontal component is retained, and the short-interval eastward displacement increment and short-interval northward displacement increment are obtained. S23, perform heading calculation and quadrant correction on the short-interval eastward displacement increment and short-interval northward displacement increment to obtain the short-interval reference heading, and store it in the short-interval buffer as a verification benchmark for heading integrity detection.

[0009] Further, in step S22, the relative position increment in the short-interval geocentric geofixed coordinate system is transformed to the ship's local northeast-sky coordinate system, separating the vertical undulation displacement and retaining the horizontal component, including: S221, obtain the ship's latitude and longitude coordinates at the current epoch, and construct the rotation matrix from the geocentric earth-fixed coordinate system to the ship's local northeast-sky coordinate system; S222, using the rotation matrix, the relative position increment in the short-interval geocentric geofixed coordinate system is transformed to the ship's local northeast-sky coordinate system to obtain the eastward displacement component, the northward displacement component, and the skyward displacement component. S223, remove the celestial displacement component, and use the eastward displacement component and the northward displacement component as the short-interval eastward displacement increment and the short-interval northward displacement increment, respectively.

[0010] Furthermore, step S3, which involves selecting historical epochs from the long-interval buffer that satisfy at least one of a preset time threshold and a preset horizontal distance threshold as reference epochs, includes: S311, traverse all historical epochs in the long interval buffer, and filter out epochs whose horizontal relative distance from the current epoch is greater than or equal to the preset horizontal distance threshold, or whose time difference from the current epoch is greater than or equal to the preset time threshold, as a set of candidate reference epochs. S312, determine whether the candidate reference epoch set is empty; if it is not empty, select the epoch with the shortest time interval with the current epoch from the candidate reference epoch set as the reference epoch. S313, if the candidate reference epoch set is empty, then the historical epoch with the largest time difference from the current epoch in the long interval buffer is selected as the reference epoch.

[0011] Further, after performing ephemeris orbital error compensation and atmospheric delay compensation on the current epoch and the reference epoch, step S3, which involves performing cross-epoch carrier phase differential calculation to obtain the long-interval coarse heading, includes: S321, Read the satellite broadcast ephemeris corresponding to the current epoch and the reference epoch respectively, calculate the satellite orbit coordinates and satellite clock difference at the two epochs, and calculate the line-of-sight orbit residual error; S322 uses the ionospheric model and the tropospheric model to calculate the atmospheric delay difference between two epochs; S323, using the line-of-sight orbital residual error and the atmospheric delay difference, pre-correct the carrier phase observation values ​​of the current epoch and the reference epoch; S324, perform cross-epoch difference on the pre-corrected carrier phase observations, solve for the long-interval relative position increment, and obtain the long-interval coarse heading through coordinate transformation and heading calculation.

[0012] Further, step S321, which involves calculating the satellite orbit coordinates and satellite clock error at two epochs and calculating the line-of-sight orbit residual error, includes: S3211, Read the satellite broadcast ephemeris corresponding to the current epoch and the reference epoch respectively, and calculate the satellite orbit coordinates and satellite clock difference at the two epochs; S3212, calculates the satellite orbital position deviation and satellite clock error difference between two epochs; S3213, Project the satellite orbital position deviation and the satellite clock error difference onto the satellite-receiver line-of-sight direction to obtain the line-of-sight orbital residual error.

[0013] Further, step S322, which involves using an ionospheric model and a tropospheric model to calculate the atmospheric delay difference between two epochs, includes: S3221, Using the Klobuchar ionospheric model, the ionospheric delay between the current epoch and the reference epoch is calculated to obtain the ionospheric delay difference. S3222, using the Saastamoinen tropospheric model, calculate the tropospheric delay between the current epoch and the reference epoch to obtain the tropospheric delay difference; S3223, The ionospheric delay difference value and the tropospheric delay difference value are combined to obtain the atmospheric delay difference value.

[0014] Further, step S324, which involves performing cross-epoch differential analysis on the pre-corrected carrier phase observations to solve for the long-interval relative position increments and obtaining the long-interval coarse heading through coordinate transformation and heading calculation, includes: S3241 performs cross-epoch difference on the pre-corrected carrier phase observations to solve for the relative position increment in the long-interval geocentric-geofixed coordinate system. S3242, the relative position increment under the long-interval geocentric geofixed coordinate system is transformed to the ship's local northeast-sky coordinate system, and the vertical undulation displacement is separated to obtain the long-interval eastward displacement increment and the long-interval northward displacement increment. S3243, the long-interval coarse heading is obtained by calculating the heading of the long-interval eastward displacement increment and the long-interval northward displacement increment using the arctangent function.

[0015] Further, step S4, which uses the short-interval reference heading as a verification benchmark to perform heading autonomy integrity testing on the long-interval coarse heading to obtain a stable heading output for the current epoch, includes: S41, extract the short-interval reference heading from the short-interval buffer, obtain the long-interval coarse heading, and calculate the absolute value of the heading difference between the two; S42, compare the absolute value of the heading difference with the preset integrity threshold; if the absolute value of the heading difference is less than or equal to the preset integrity threshold, then take the long-interval coarse heading as the effective output heading; S43, if the absolute value of the heading difference is greater than the preset integrity threshold, it is determined that the long-interval coarse heading is abnormal, and the short-interval reference heading is output as the alternative heading for the current epoch. S44. Perform quadrant correction and normalization on the output heading determined in step S42 or step S43 to obtain stable heading data for the current epoch.

[0016] Furthermore, after obtaining the stable heading output for the current epoch in step S4, the process also includes: S51, obtain the speed and heading values ​​for a consecutive preset number of observation epochs, and calculate the speed fluctuation value and heading change rate; S52, Based on the comparison results of the speed fluctuation value and the heading change rate with the preset speed fluctuation threshold and the preset heading change rate threshold, determine the current navigation condition of the ship. S53 adaptively adjusts the differential strategy and detection threshold based on the determined navigation conditions.

[0017] Furthermore, step S53, which involves adaptively adjusting the differential strategy and detection threshold based on the determined navigation conditions, includes: S531, if it is determined that the ship is in a uniform speed straight navigation condition, continue to execute steps S2 to S4 in a loop, and use the long-interval coarse heading as the main output heading. S532, if it is determined that the ship has entered a turning or speed change maneuvering condition, the preset time threshold, the preset horizontal distance threshold and the preset integrity threshold are increased, and the maximum time difference of the reference epochs that can be screened in the long interval buffer is shortened; S533, under the steering or speed change maneuvering condition, reduce the output weight of the long-interval coarse heading and increase the output weight of the short-interval reference heading.

[0018] Accordingly, a second aspect of the present invention provides an electronic device, comprising: at least one processor; and a memory connected to the at least one processor; wherein the memory stores instructions executable by the at least one processor, the instructions being executed by the at least one processor to cause the at least one processor to perform the above-described method for measuring the anti-turbulence heading stability of marine GNSS based on cross-epoch difference.

[0019] Accordingly, a third aspect of the present invention provides a computer-readable storage medium having computer instructions stored thereon, which, when executed by a processor, implement the above-described method for measuring the anti-turbulence heading stability of marine GNSS based on cross-epoch difference.

[0020] The above-described technical solutions of the embodiments of the present invention have the following beneficial technical effects: 1. By constructing a two-layer architecture that coordinates short-interval differential and long-interval cross-epoch differential operations, and using the short-interval differential heading as the verification benchmark, this system innovatively achieves receiver autonomous integrity monitoring in the heading dimension. The short-interval differential automatically eliminates integer ambiguity by utilizing the carrier phase difference between adjacent epochs, and can output a reference heading without accumulated error in real time without relying on external sensors. The long-interval differential, by selecting historical epochs that meet time or distance thresholds, uses the long baseline relative displacement vector to characterize the ship's macroscopic navigation direction, which significantly dilutes the instantaneous position and velocity disturbance errors introduced by sea surface wind and wave turbulence from a mechanistic perspective. On this basis, by comparing the consistency between the short-interval reference heading and the long-interval coarse heading, when the long-interval result is abnormal due to cycle slip, multipath, or residual orbit error, it can autonomously determine and switch to the short-interval heading output, fundamentally solving the defect of traditional single-antenna direction finding schemes lacking integrity monitoring capabilities at the heading level, and significantly improving the reliability of heading output. 2. To address the technical challenges of systematic displacement deviations caused by the accumulation of satellite broadcast ephemeris orbital bias, satellite clock bias, and ionospheric and tropospheric delays over long time in long-time-difference cross-epoch differential scenarios, this paper proposes a method that proactively reads the broadcast ephemeris corresponding to the current epoch and the reference epoch before performing long-interval cross-epoch differential calculations. This method calculates the satellite orbital coordinates and satellite clock bias at both times, calculates the line-of-sight orbital residual error, and simultaneously compensates for the atmospheric delay difference between the two epochs using ionospheric and tropospheric models. After pre-correcting the carrier phase observations using all the aforementioned error terms, the differential calculation is performed. This compensation process effectively eliminates the impact of long-timescale systematic displacement shifts on relative position increments, improves the upper limit of the calculation accuracy of long-interval differentials, and enables high-precision heading measurement results even with significantly extended time intervals, thus broadening the applicable time range for cross-epoch differential calculations. 3. By monitoring the ship's navigation status in real time, the system automatically identifies uniform straight-line navigation and turning / variable-speed maneuvering conditions, and adaptively adjusts the differential strategy and detection thresholds according to the condition. During uniform cruising or low-speed berthing, long-interval differential results are prioritized as the primary output heading, fully leveraging their advantages in anti-turbulence and high stability. During turning or variable-speed maneuvering, the system automatically increases the time threshold, distance threshold, and integrity detection threshold, shortens the allowable time difference upper limit of long-interval reference epochs, reduces the weight of long-interval differential output, and increases the weight of short-interval differential output, thus ensuring both heading response speed and output stability. This adaptive mechanism resolves the inherent contradiction that fixed differential strategies cannot simultaneously satisfy high heading stability and high real-time performance, enabling this solution to adapt to the heading measurement needs of ships in various sea states and all operating conditions, possessing universal capability across all scenarios. Attached Figure Description

[0021] Figure 1This is a flowchart of a method for measuring the heading stability of a shipborne GNSS based on cross-epoch difference, provided in an embodiment of the present invention. Detailed Implementation

[0022] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to specific embodiments and the accompanying drawings. It should be understood that these descriptions are merely exemplary and not intended to limit the scope of the invention. Furthermore, descriptions of well-known structures and techniques are omitted in the following description to avoid unnecessarily obscuring the concept of the invention.

[0023] During maritime navigation, especially in low-speed maneuvering scenarios such as berthing in port, fixed-point operations, and channel fine-tuning, ships require real-time, stable, and reliable heading information to support control decisions. However, the hull rolling caused by wind and waves leads to severe headway fluctuations in traditional direction-finding methods based on single-antenna GNSS receivers, seriously affecting navigation safety and operational accuracy. Existing solutions, such as heading calculation methods based on Doppler instantaneous velocity, cannot distinguish between the ship's actual movement and the speed of wind and wave disturbances, resulting in heading fluctuations of several degrees or even more than ten degrees at low speeds. While orientation methods based on dual-antenna baselines can provide relatively stable headings, the hardware cost increases exponentially and installation requirements are stringent, making them difficult to implement on small and medium-sized vessels. Therefore, achieving highly stable heading measurement with resistance to turbulence under a single-antenna hardware architecture has become a pressing technical problem to be solved in this field.

[0024] To address the aforementioned issues, this embodiment proposes a method for measuring the turbulence-resistant heading stability of marine GNSS based on cross-epoch differential methods, constructing a two-layer measurement architecture that coordinates short-interval differential and long-interval differential methods. This architecture establishes two independent sliding buffer queues at the data layer, providing adjacent epoch data for the short-interval differential method and historical epoch data for the long-interval differential method. At the solution layer, the short-interval differential method quickly solves the instantaneous reference heading through adjacent epoch carrier phase differential, while the long-interval differential method, after selecting suitable distant reference epochs and compensating for orbital errors and atmospheric delays, solves for the long baseline displacement vector heading with natural turbulence resistance. At the output layer, using the short-interval reference heading as a verification benchmark, autonomous integrity checks are performed on the heading dimension of the long-interval coarse heading, and automatic identification of abnormal headings and optimal output of reliable headings are achieved through consistency comparison. This architecture eliminates the need for additional hardware such as inertial measurement units, electronic compasses, or dual antennas. It relies entirely on the raw observation data from a single-antenna GNSS receiver to complete all calculations, simultaneously improving the stability and reliability of the heading from both the measurement mechanism and data verification perspectives.

[0025] Please refer to Figure 1The first aspect of this invention provides a method for measuring the anti-turbulence heading stability of marine GNSS based on cross-epoch difference, comprising the following steps: S1: Real-time acquisition of raw observation data from the shipborne GNSS receiver; construction of a two-layer sliding buffer queue; the two-layer sliding buffer queue includes a long-interval buffer for caching observation data from multiple historical epochs, and a short-interval buffer for caching observation data from the current epoch and the previous adjacent epoch.

[0026] The shipborne GNSS receiver continuously outputs raw observation data for each observation epoch at a preset sampling frequency. This raw observation data includes at least the epoch timestamp, satellite pseudorange observations, carrier phase observations, carrier-to-noise ratio information, and latitude and longitude coordinates calculated by the receiver. The receiver is typically installed in an open location such as the ship's superstructure or mast to ensure continuous and stable reception of satellite signals. To meet the parallel processing requirements of subsequent long and short epoch differential processing, two independent sliding buffer queues are constructed at the receiver's data output. The long-interval buffer stores complete observation data from multiple consecutive historical epochs. Its buffer duration can be set according to the ship's typical sailing speed and sea conditions. For example, in applications on small and medium-sized ships, a buffer depth of no less than 30 seconds can be set to ensure that reference epochs meeting the spatial interval requirements can still be selected even when the ship is hovering at low speed. The short-interval buffer retains only the observation data from the current epoch and the immediately preceding epoch, used for real-time short-timescale differential preprocessing. Both cache queues employ a first-in, first-out (FIFO) sliding update mechanism. Whenever new observation history data arrives, it is synchronously written to both caches, while historical data exceeding the cache duration is moved out of the long-interval cache. This two-layer cache architecture provides an independent and non-interfering data foundation for the subsequent simultaneous execution of short-interval differential verification and long-interval differential calculation.

[0027] In a preferred embodiment of this example, the shipborne GNSS receiver continuously acquires raw observation data at a frequency of 10Hz. The buffer duration of the long-interval buffer is set to be no less than 30 seconds to ensure that historical reference epochs meeting the spatial interval requirements can still be selected even when the ship is lingering at low speed. In step S3, the preset horizontal distance threshold is set to 10 meters, and the preset time threshold is set to 5 seconds. Meeting either condition will allow the corresponding historical epoch to be included in the candidate reference epoch set. In step S4, the base value of the preset integrity threshold is set to 2 degrees, which can be adaptively adjusted to 3 to 5 degrees under high sea state conditions to accommodate the reasonable deviation range caused by the increased fluctuation of the short-interval reference heading under severe sea states. In step S5, the preset speed fluctuation threshold is set to 1.0 knots, the preset heading change rate threshold is set to 0.5 degrees per second, and the number of observation epochs required for continuous determination is set to 3.

[0028] S2, perform carrier phase difference between the current epoch and the previous adjacent epoch in the short interval buffer to calculate the short interval reference heading.

[0029] Using the carrier phase observations of the current epoch and the previous adjacent epoch stored in the short-interval buffer, carrier phase differential calculations are performed between adjacent epochs. Since the time interval between adjacent epochs is extremely short, typically on the order of the receiver sampling period, the integer ambiguity parameter of the carrier phase remains unchanged within this timescale. Integer ambiguity terms can be automatically eliminated by subtracting between epochs, without the need for complex ambiguity fixation calculations. The short baseline relative position increment vector in the geocentric-fixed coordinate system is obtained after differential calculation. Subsequently, a rotation matrix is ​​constructed using the ship's latitude and longitude coordinates of the current epoch to transform this relative position increment from the geocentric-fixed coordinate system to a local northeast-sky coordinate system centered on the ship. In the transformed northeast-sky component, the celestial displacement component mainly reflects the vertical undulation motion of the ship under the action of waves and is not directly related to the ship's horizontal heading; therefore, it is separated and removed, retaining only the eastward and northward displacement components. Based on the aforementioned horizontal displacement components, a short-interval reference heading angle is calculated using arctangent calculation. Quadrant correction is then performed by combining the positive and negative signs of the eastward and northward displacement components, normalizing the heading angle to the range of 0° to 360°. This short-interval reference heading is stored in the short-interval buffer as a verification benchmark for subsequent heading integrity testing. Its advantages lie in its extremely short time baseline, the absence of error accumulation effects, and its ability to accurately reflect the direction of motion measured by the receiver at the current instant.

[0030] S3: Select historical epochs from the long-interval buffer that meet at least one of the preset time threshold and preset horizontal distance threshold as reference epochs. After performing ephemeris orbital error compensation and atmospheric delay compensation on the current epoch and the reference epoch, perform cross-epoch carrier phase differential to calculate the long-interval coarse heading.

[0031] Suitable reference epochs are selected from multiple historical epochs stored in a long-interval buffer. The selection criteria employ a parallel approach using time thresholds and horizontal distance thresholds. Specifically, each historical epoch in the long-interval buffer is traversed, and its time difference with the current epoch and the horizontal relative distance between the receiver positions of the two epochs are calculated. If a historical epoch satisfies either a time difference greater than or equal to a preset time threshold, or a horizontal relative distance greater than or equal to a preset horizontal distance threshold, it is included in the candidate reference epoch set. When multiple candidate epochs exist, the epoch with the shortest time interval to the current epoch is selected as the differential reference, minimizing the time span while satisfying the spatial separation condition. This selection mechanism considers scenarios that may occur during actual ship navigation, such as prolonged periods of stillness or extremely low-speed hovering. By using the time threshold as a supplementary condition to the spatial threshold, it ensures that effective long-interval reference epochs can be selected under various operating conditions. After selecting the reference epochs, carrier phase differential is not performed directly; instead, a system error compensation step is executed first. Specifically, the satellite broadcast ephemeris parameters corresponding to the current epoch and the reference epoch are read respectively, and the satellite's three-dimensional orbital coordinates and satellite clock errors at the two epochs are calculated. The satellite orbital position deviation and satellite clock error difference between the two epochs are calculated, and the above deviations are projected onto the line-of-sight direction between the satellite and the receiver to obtain the line-of-sight orbital residual error. At the same time, the atmospheric propagation delay of the two epochs is calculated using the ionospheric delay model and the tropospheric delay model respectively, to obtain the ionospheric delay difference and the tropospheric delay difference. Using the line-of-sight orbital residual error and the atmospheric delay difference, the carrier phase observations of the two epochs are pre-corrected to eliminate the systematic deviations introduced by the accumulation of orbital errors and atmospheric delays under long timescale intervals. After compensation, cross-epoch difference is performed on the corrected carrier phase observations to solve for the relative position increment in the long-interval geocentric Earth-fixed coordinate system. After coordinate transformation to the local Northeast-Sky coordinate system, the long-interval eastward displacement increment and northward displacement increment are obtained, and then the long-interval coarse heading is calculated. Because the reference epoch selected by the long-interval differential has a large time span or spatial distance between it and the current epoch, its relative displacement vector is mainly contributed by the actual navigation motion of the ship. The instantaneous disturbance caused by sea surface turbulence is greatly diluted in the long baseline. Therefore, the long-interval coarse heading has natural anti-turbulence stability characteristics.

[0032] When selecting reference epochs from the long-interval buffer, this embodiment employs a three-tiered progressive selection logic: "parallel selection with dual thresholds, optimal selection from multiple candidates, and fallback guarantee without candidates." First, all historical epochs stored in the long-interval buffer are traversed, and the time difference and horizontal relative distance between each historical epoch and the current epoch are calculated. Historical epochs that satisfy either a time difference greater than or equal to a preset time threshold or a horizontal relative distance greater than or equal to a preset horizontal distance threshold are included in the candidate reference epoch set. Second, it is determined whether the candidate reference epoch set is empty. If the set is not empty, it indicates that at least one candidate epoch satisfies the dual threshold conditions. In this case, the epoch with the shortest time interval to the current epoch is selected as the final long-interval differential reference epoch. This minimizes the time span while satisfying the spatial separation condition, thereby reducing the residual amount of orbital error and atmospheric delay compensation. If the candidate reference epoch set is empty, for example, if the buffer data is not yet full during the initialization phase after the receiver starts up, or if the ship is stationary for a long time and the horizontal distance cannot reach the preset threshold, then the historical epoch with the largest time difference from the current epoch in the long interval buffer is directly selected as the reference epoch, that is, the earliest epoch stored in the buffer is selected to ensure the continuity of the long interval differential process.

[0033] S4, using the short-interval reference heading as the verification benchmark, performs autonomous heading integrity testing on the long-interval coarse heading to obtain a stable heading output for the current epoch. The autonomous heading integrity testing includes: calculating the heading difference between the short-interval reference heading and the long-interval coarse heading; if the heading difference exceeds a preset integrity threshold, the long-interval coarse heading is determined to be abnormal and the short-interval reference heading is output; otherwise, the long-interval coarse heading is output.

[0034] Using the short-interval reference heading obtained in step S2 as a redundancy check benchmark, an autonomous integrity check is performed on the long-interval coarse heading obtained in step S3. First, the short-interval reference heading and the long-interval coarse heading corresponding to the current epoch are extracted, and the absolute value of the heading difference between them is calculated. This heading difference is compared with a preset heading integrity threshold. If the heading difference is less than or equal to the preset integrity threshold, it indicates that the long-interval differencing result has good consistency with the short-interval reference heading, and the long-interval differencing process is not significantly interfered with by abnormal factors such as cycle slips, multipath effects, or residual orbital errors. The long-interval coarse heading is then determined to have passed the integrity check and is used as the valid output heading for the current epoch. If the heading difference exceeds the preset integrity threshold, it indicates that the long-interval differencing result is abnormal. Possible causes include, but are not limited to, cycle slips in the carrier phase of the reference epoch or the current epoch, residual errors exceeding the expected range during the long-interval differencing process, or the observation value of a certain satellite being severely affected by multipath effects. At this point, the heading integrity alarm mechanism is triggered, abandoning the long-interval coarse heading and instead outputting a short-interval reference heading as the alternative heading for the current epoch, thereby preventing erroneous heading information from being output to the ship's control system. After determining the output heading, quadrant correction and normalization are performed, uniformly standardizing the heading angle to a 0° to 360° coordinate system with true north as the 0° reference, as the final stable heading data output for the current epoch.

[0035] Through the above steps, a two-layer measurement architecture was constructed using only a single-antenna GNSS receiver, with short-interval differential and long-interval differential working in tandem. The short-interval differential provides a reference heading with high real-time performance and no accumulated error, while the long-interval differential utilizes the long baseline displacement vector to dilute sea surface turbulence disturbances, resulting in a highly stable coarse heading. By performing autonomous integrity detection of the long-interval coarse heading using the short-interval reference heading, anomaly identification and automatic correction are completed before heading output. This achieves a simultaneous improvement in heading measurement stability and output reliability under low-speed turbulence conditions, effectively solving the technical problems of severe heading jitter and lack of integrity monitoring methods in traditional single-antenna direction finding schemes.

[0036] Specifically, step S2 calculates the short-interval reference heading, including: S21, perform epoch-to-epoch difference on the carrier phase observations of the current epoch and the previous adjacent epoch to eliminate carrier integer ambiguity and solve for the relative position increment in the short-interval geocentric geofixed coordinate system.

[0037] Carrier phase observations include the geometric distance from the satellite to the receiver, receiver clock bias, satellite clock bias, atmospheric delay, and integer ambiguity. The integer ambiguity remains constant under conditions of continuous satellite signal tracking and no cycle slips. By subtracting the carrier phase observations from adjacent epochs, all error terms that remain unchanged for a very short time, including integer ambiguity, are automatically eliminated, eliminating the need for complex ambiguity searches and fixed calculations, greatly simplifying the process of calculating relative position increments. In practice, common-view satellites at two epochs can be matched first, eliminating satellites with observations only in a single epoch. For each commonly observed satellite, epoch-level carrier phase difference observations are constructed separately. Then, the least squares method or weighted least squares method is used to solve for the three-dimensional relative position increment vector in the geocentric-ground-fixed coordinate system. Since the time interval between adjacent epochs is usually only the receiver sampling interval, the carrier phase observation noise remains low after epoch-level differencing, resulting in high accuracy for the short-interval relative position increments, providing a reliable data foundation for subsequent heading calculations.

[0038] S22 transforms the relative position increment in the short-interval geocentric-ground-fixed coordinate system to the ship's local northeast-sky coordinate system, separates the vertical undulation displacement, retains the horizontal component, and obtains the short-interval eastward displacement increment and the short-interval northward displacement increment.

[0039] The geocentric-fixed coordinate system is a three-dimensional rectangular coordinate system with the Earth's center of mass as its origin. Its coordinate axes are fixed to the Earth. While convenient for unified description of satellite orbits, it does not directly reflect the ship's direction of motion in the horizontal plane. Therefore, it is necessary to transform the relative position increment to a local northeast-sky coordinate system with the ship's current position as its origin. This coordinate system uses the east direction as the positive X-axis, the north direction as the positive Y-axis, and the sky direction as the positive Z-axis. The heading angle is defined as the angle rotated clockwise from the north direction to the ship's horizontal direction of motion. The transformation process first uses the latitude and longitude coordinates obtained from the current epoch receiver to calculate the rotation matrix from the geocentric-fixed coordinate system to the local northeast-sky coordinate system. Then, the relative position increment vector is multiplied by this rotation matrix to obtain three mutually orthogonal displacement components: east, north, and sky. Based on this, the sky displacement component is discarded, retaining only the east and north horizontal displacement components. The vertical displacement component is mainly contributed by factors such as the heave motion of the hull under the action of waves and the changes in tidal level. It is unrelated to the ship's heading in the horizontal plane. Separating and eliminating it can effectively suppress the interference of vertical turbulence motion on the heading calculation, so that the heading angle is determined only by the actual motion displacement in the horizontal plane.

[0040] S23, calculate the heading and quadrant correction for the short-interval eastward displacement increment and the short-interval northward displacement increment to obtain the short-interval reference heading, and store it in the short-interval buffer as the verification benchmark for heading integrity detection.

[0041] The basic formula for calculating the heading angle is the arctangent function of the ratio of the northward displacement increment to the eastward displacement increment. This calculation yields an angle ranging from -90 degrees to +90 degrees, covering only two quadrants. To obtain a heading angle continuously expressed from 0 degrees to 360 degrees relative to true north, quadrant correction is required based on the signs of the eastward and northward displacement increments. Specifically, when the northward displacement increment is positive, the heading angle is in the first or second quadrant; when the northward displacement increment is negative, the heading angle is in the third or fourth quadrant. Furthermore, the signs of the eastward displacement increment are used to further distinguish between the first and second quadrants, and between the third and fourth quadrants. The quadrant-corrected heading angle fully covers the range from 0 degrees to 360 degrees, accurately representing the ship's horizontal movement direction relative to true north. The heading angle is stored in the short-interval buffer as a short-interval reference heading. The real-time heading estimate obtained based on the shortest time baseline has almost no error accumulation effect due to the extremely small epoch interval. It can reflect the horizontal motion direction perceived by the receiver at the current moment in a relatively realistic way. Therefore, it can be used as a reliable benchmark to verify whether the long-interval differential result has abnormal drift in subsequent steps.

[0042] Through the above steps, the short-interval differential process utilizes the characteristic of automatically eliminating integer ambiguity by using phase difference between adjacent epoch carriers, avoiding the complex ambiguity search and fixing process in traditional positioning solutions, and obtaining high-precision relative position increments with minimal computational overhead. Subsequently, through coordinate transformation and vertical separation, the heave motion component of the hull under the action of waves is separated from the heading calculation link, so that the heading angle is determined only by the actual displacement in the horizontal plane, suppressing the interference of vertical turbulence on heading measurement from a mechanistic perspective. The final output short-interval reference heading has the characteristics of no cumulative error and high real-time performance, providing a reliable verification benchmark independent of long-interval differential for subsequent heading integrity detection, enabling the entire system to autonomously complete the identification and judgment of heading anomalies even in the absence of external reference information.

[0043] Further, in step S22, the relative position increments in the short-interval geocentric geofixed coordinate system are transformed to the ship's local northeast-sky coordinate system, separating the vertical undulation displacement while retaining the horizontal component, including: S221, obtain the ship's latitude and longitude coordinates at the current epoch, and construct the rotation matrix from the geocentric geofixed coordinate system to the ship's local northeast-sky coordinate system.

[0044] The geocentric-fixed coordinate system is a three-dimensional rectangular coordinate system with the Earth's center of mass as its origin and fixed to the Earth. Its coordinate axes are defined by the Earth's rotation axis and the prime meridian. Satellite orbit coordinates and receiver positions are typically expressed in this coordinate system. However, ship heading, as a horizontal directional quantity, only has intuitive physical meaning in a local horizontal coordinate system centered on the ship's current position. Therefore, converting the relative position increment in the geocentric-fixed coordinate system into displacement components in the local northeast-sky coordinate system is a necessary prerequisite for subsequent heading calculations. The process of constructing the rotation matrix relies on the ship's geodetic coordinates at the current epoch. Specifically, using the longitude and latitude of the ship's location, the direction cosines of the three coordinate axes of the geocentric-fixed coordinate system in the local northeast-sky coordinate system are determined, thus forming a three-row, three-column orthogonal rotation matrix. Each row of this rotation matrix corresponds to the unit vector expression of the eastward, northward, and celestial axes of the local northeast-sky coordinate system in the geocentric-fixed coordinate system. In actual calculations, the elements of the rotation matrix are only related to the latitude and longitude of the ship at the current epoch. Therefore, the rotation matrix needs to be reconstructed based on the real-time position for each observation epoch to adapt to the position changes of the ship in different navigation areas and at different times.

[0045] S222, using a rotation matrix, transforms the relative position increment in the short-interval geocentric geofixed coordinate system to the ship's local northeast-sky coordinate system, obtaining the eastward displacement component, the northward displacement component, and the skyward displacement component.

[0046] The three-dimensional relative position increment vector in the short-interval geocentric-fixed coordinate system obtained in step S21 is multiplied by the rotation matrix constructed in step S221 to complete the linear transformation from the geocentric-fixed coordinate system to the ship's local northeast-sky coordinate system. The displacement vector originally expressed in the coordinate system with the Earth's center of mass as the origin is reprojected into the local horizontal coordinate system with the ship's current position as the origin, thereby decomposing the ship's motion components in the horizontal and vertical directions. Of the three displacement components obtained after the transformation, the eastward displacement component represents the ship's position change in the east-west direction, the northward displacement component represents the ship's position change in the north-south direction, and the skyward displacement component represents the ship's position change in the vertical direction. In typical ship navigation scenarios at sea, the eastward and northward displacement components are mainly contributed by factors such as the ship's horizontal navigation motion and wind and pressure differences, while the skyward displacement component mainly comes from the ship's heave motion caused by ocean waves and water level changes caused by tides. The separation of these three components provides a data basis for the subsequent selective retention of the horizontal component and the elimination of vertical interference.

[0047] S223, remove the celestial displacement component, and take the eastward displacement component and the northward displacement component as the short-interval eastward displacement increment and the short-interval northward displacement increment, respectively.

[0048] The three-dimensional displacement components obtained after coordinate transformation in step S222 are selectively processed. The upward displacement component is discarded, and only the eastward and northward displacement components are retained as effective horizontal displacement increments. The ship's heading is defined as the angle in the horizontal plane rotated clockwise from due north to the ship's horizontal direction of motion. Its calculation depends only on the eastward and northward horizontal displacement components; the upward displacement component does not participate in the geometric calculation of the heading angle. More importantly, under rough sea conditions, the heave motion of the hull under the action of waves can reach several meters. There is no definite geometric relationship between this vertical motion and the ship's horizontal heading. If the upward component is not separated and discarded, the vertical displacement will indirectly affect the calculation accuracy of the horizontal component through the coupling effect of coordinate transformation, thus introducing additional heading errors. Furthermore, in short-interval difference scenarios, due to the extremely short time interval between adjacent epochs, the horizontal displacement increment itself is small, and the proportion of vertical disturbance in the relative displacement is relatively large, making the separation and discarding effect particularly significant. Therefore, by eliminating the axial displacement component, the interference of the ship's heave motion on the heading calculation link can be effectively eliminated. This ensures that the short-interval reference heading obtained in subsequent calculations only reflects the ship's true direction of motion in the horizontal plane, significantly improving the accuracy and stability of the short-interval reference heading under turbulent sea conditions. The eastward and northward displacement increments after the above processing are passed to subsequent heading calculation steps, completing the dimensionality reduction from three-dimensional displacement to two-dimensional horizontal displacement, providing clean input data for accurate heading angle calculation.

[0049] Specifically, step S3, which involves selecting historical epochs from the long-interval buffer that satisfy at least one of a preset time threshold and a preset horizontal distance threshold as reference epochs, includes: S311, traverse all historical epochs in the long-interval buffer, and select epochs whose horizontal relative distance from the current epoch is greater than or equal to a preset horizontal distance threshold, or whose time difference from the current epoch is greater than or equal to a preset time threshold, as candidate reference epoch sets.

[0050] The entire historical observation epoch stored in the long-interval buffer is iterated through one by one. The time interval and horizontal relative distance between each historical epoch and the current epoch are calculated. Candidate epochs are initially screened based on preset dual threshold conditions. The time interval refers to the absolute value of the difference between the timestamp of the current epoch and the timestamp of a historical epoch. The horizontal relative distance refers to the Euclidean distance in the horizontal plane calculated using the latitude and longitude coordinates of the receivers corresponding to the two epochs. The preset time threshold and preset horizontal distance threshold are determined in parallel; that is, a historical epoch can be included in the candidate reference epoch set if it meets either of the two conditions. This parallel screening mechanism considers the actual needs of the ship under different navigation states: when the ship is traveling at normal speed, sufficient horizontal displacement can be accumulated in a short time, and the horizontal distance threshold condition will be met first; however, when the ship is in a low-speed loitering or stationary operation state, the horizontal displacement accumulates slowly, and relying solely on the distance threshold may result in a long time without selecting a reference epoch. In this case, the time threshold is introduced as a supplementary condition; even if the horizontal distance has not yet reached the preset requirement, as long as the time interval is long enough, the historical epoch can still be considered as a candidate. The specific values ​​of the two thresholds can be set according to the typical speed of the ship and sea conditions. For example, for small and medium-sized ships with low speeds, the horizontal distance threshold can be set to 10 meters and the time threshold can be set to 5 seconds to ensure that effective candidate reference eras can be selected in a timely manner under various working conditions.

[0051] S312, determine whether the candidate reference epoch set is empty; if it is not empty, select the epoch with the shortest time interval with the current epoch from the candidate reference epoch set as the reference epoch.

[0052] The candidate reference epochs obtained in step S311 are evaluated for existence and then selected based on their merit. First, it is determined whether the candidate reference epoch set is empty. If the set is not empty, it indicates that at least one historical epoch in the long-interval buffer meets the double-threshold screening condition. Then, among the multiple candidate epochs, the epoch with the shortest time interval to the current epoch is selected as the final long-interval differential reference epoch. Under the premise of satisfying the spatial separation condition, the shorter the time interval, the smaller the accumulated amount of slowly varying errors such as orbital errors of satellite broadcast ephemeris, satellite clock errors, and atmospheric delays between the two epochs. This also results in smaller residual errors during subsequent systematic error compensation, which is beneficial for obtaining higher-precision long-interval relative displacement calculation results. Simultaneously, a shorter time interval also means that the ship's trajectory between two observations is simpler, the possibility of significant changes in heading is lower, and the calculated heading is more representative of the ship's actual navigation direction at the current moment. In practice, historical epochs in the long-interval cache can be sorted by timestamp, and the search can be performed backward from the candidate epoch closest to the current epoch. The first epoch that meets the double threshold condition is the candidate epoch with the shortest time interval, thereby avoiding repeated traversal of all cached epochs and improving the filtering efficiency.

[0053] S313 If the candidate reference epoch set is empty, then the historical epoch with the largest time difference from the current epoch in the long interval buffer is selected as the reference epoch.

[0054] An empty candidate reference epoch set typically occurs when a ship is stationary or drifting at extremely low speeds for an extended period. This occurs in an extreme combination where, within the entire buffer duration of the long-interval buffer, no historical epoch meets the preset horizontal distance threshold, and the time difference between all historical epochs also fails to reach the preset time threshold. Alternatively, it may occur during the initialization phase when the buffer is not yet full, and the total amount of cached historical data is insufficient to cover the time span required by the preset time threshold. In such cases, strictly adhering to the dual-threshold criteria will fail to select any reference epochs, leading to the interruption of the long-interval difference process. To ensure the continuity and robustness of the heading measurement method, this invention employs a fallback selection strategy: directly selecting the historical epoch with the largest time difference from the current epoch in the long-interval buffer as the reference epoch, i.e., selecting the earliest historical epoch stored in the buffer. When spatial separation conditions cannot be met, selecting the epoch with the longest time interval maximizes the positional difference between the two epochs. Even if this difference is not caused by the ship's active navigation, the accumulated small drift can form a usable baseline length over a long timescale, thus ensuring the basic usability of the long-interval difference. Meanwhile, this fallback strategy also ensures that the method can output heading results normally during the initial buffer setup phase of receiver startup, without having to wait for the long interval buffer to be completely filled before starting long interval differential operations.

[0055] Through the above steps, the screening process of long-interval reference epochs forms a three-level progressive selection logic of "parallel screening with dual thresholds, selection of the best candidate from multiple candidates, and no candidate fallback guarantee". It can stably select suitable long-interval reference epochs under various working conditions such as normal ship navigation, low-speed wandering, stationary berthing, and equipment initialization, providing a reliable epoch benchmark for subsequent high-precision solution of long-interval difference.

[0056] Specifically, in step S3, after performing ephemeris orbital error compensation and atmospheric delay compensation between the current epoch and the reference epoch, cross-epoch carrier phase differential is performed to calculate the long-interval coarse heading, including: S321 reads the satellite broadcast ephemeris corresponding to the current epoch and the reference epoch respectively, calculates the satellite orbit coordinates and satellite clock difference at the two epochs, and calculates the line-of-sight orbit residual error.

[0057] Broadcast ephemeris is a set of orbital prediction parameters broadcast in real time by GNSS satellites. Receivers can obtain the satellite's orbital position and clock bias information at a specified time by decoding the ephemeris message. Because the orbital prediction accuracy of broadcast ephemeris is limited, typically on the order of meters to several meters, and the prediction error gradually increases over time, the difference in orbital position deviation and satellite clock bias between two epochs becomes significant when the time interval between two epochs is long. First, the three-dimensional orbital coordinates of each common-view satellite at the current epoch and the reference epoch are calculated using the broadcast ephemeris parameters, while simultaneously extracting the satellite clock bias parameters at both epochs. Then, for each common-view satellite, the satellite orbital position deviation vector and the satellite clock bias difference between the two epochs are calculated. The orbital position deviation vector is directly taken as the difference between the satellite orbital coordinates at the two epochs, and the clock bias difference is taken as the difference between the satellite clock bias parameters at the two epochs. Based on this, the orbital position deviation vector and the clock bias difference are projected onto the line-of-sight direction between the corresponding satellite and the receiver. The line-of-sight unit vector is determined by the receiver's approximate position and the satellite's orbital coordinates. The projection operation transforms the orbital deviation in three-dimensional space into a one-dimensional distance error along the signal propagation direction, thus obtaining the line-of-sight orbital residual error. This error term directly reflects the magnitude of the impact of the broadcast ephemeris inconsistency between two epochs on the carrier phase observations and is an important input for subsequent pre-correction stages.

[0058] S322 uses an ionospheric model and a tropospheric model to calculate the atmospheric delay difference between two epochs.

[0059] When GNSS signals traverse the Earth's atmosphere, free electrons in the ionosphere and neutral gas molecules in the troposphere cause changes in signal propagation speed and path curvature, respectively, resulting in ionospheric and tropospheric delays. In short-interval differential scenarios between adjacent epochs, the atmospheric delay variation between two epochs is typically below the millimeter level and can be naturally eliminated during differential calculations due to the extremely short time interval. However, in long-interval, cross-epoch differential scenarios, the time interval between two epochs can reach several seconds or even tens of seconds. During this period, changes in satellite elevation angle, the spatiotemporal evolution of the total ionospheric electron quantity, and changes in tropospheric temperature, humidity, and pressure conditions all lead to significant changes in atmospheric delay. If not compensated, these changes will directly translate into systematic biases in the relative displacement calculation. Therefore, empirical ionospheric and tropospheric models are used to estimate the atmospheric delay at two different times. The ionospheric delay can be calculated using the Klobuchar model, which estimates the vertical ionospheric delay along the signal propagation path based on ionospheric parameters broadcast in the ephemeris, satellite elevation angle, and user location. The tropospheric delay can be calculated using the Saastamoinen model, which estimates the tropospheric zenith delay based on meteorological parameters such as surface temperature, air pressure, and water vapor partial pressure, and then projects it onto the signal propagation direction using a mapping function. After calculating the ionospheric and tropospheric delay values ​​for each common-view satellite at the current and reference epochs, the delay values ​​at the two epochs are subtracted to obtain the ionospheric and tropospheric delay differences for each satellite between the two epochs. These differences are then combined to form the atmospheric delay difference, which is used for pre-correction of subsequent carrier phase observations.

[0060] S323 uses line-of-sight orbital residual error and atmospheric delay difference to pre-correct carrier phase observations between the current epoch and the reference epoch.

[0061] The line-of-sight orbital residual error calculated in step S321 and the atmospheric delay difference calculated in step S322 are used as correction values ​​to perform satellite-by-satellite pre-correction processing on the carrier phase observations participating in long-interval cross-epoch differential analysis. The object of pre-correction can be the carrier phase observation of the current epoch, the carrier phase observation of the reference epoch, or the observations of two epochs can be allocated and corrected simultaneously. The purpose is to pre-deduct the systematic error term caused by the long time-scale interval before constructing the inter-epoch differential observation equation. Specifically, for each common-view satellite with valid observations in both epochs, its corresponding line-of-sight orbital residual error is converted into a range correction for the carrier phase. The ionospheric delay difference and the tropospheric delay difference are also converted into range corrections after being combined. Then, the above corrections are uniformly converted to the number of weeks of the carrier phase and deducted from the original carrier phase observation of the corresponding epoch. The carrier phase observations obtained after pre-correction have eliminated most of the systematic errors introduced by the inconsistency in broadcast ephemeris orbital deviations between two epochs and atmospheric delay variations. The residuals mainly originate from carrier phase observation noise, multipath effects, and the model's own approximation errors. This pre-correction step is the core element that distinguishes long-interval differencing from short-interval differencing, and it is also a key technical means to ensure that long-timescale cross-epoch differencing can still obtain high-precision relative displacement calculation results.

[0062] S324 performs cross-epoch difference on the pre-corrected carrier phase observations to solve for the long-interval relative position increment, and obtains the long-interval coarse heading through coordinate transformation and heading calculation.

[0063] First, the pre-corrected carrier phase observations at the current epoch and the reference epoch are subtracted between epochs to construct a cross-epoch carrier phase differential observation equation. Since the pre-correction step has eliminated systematic deviations caused by orbital errors and atmospheric delay variations, the main unknown parameter in the differential observation equation is the three-dimensional relative position increment of the receiver between the two epochs. Simultaneously, because the time interval between two epochs is much longer than that between adjacent epochs, the probability of cycle slips in the carrier phase over long time intervals increases significantly. In practical processing, the continuity check of the carrier phase observations and cycle slip detection markers can be combined to deweight or remove satellites experiencing cycle slips, prioritizing satellite observations that maintain continuous tracking between two epochs and have a high carrier-to-noise ratio for the solution. Subsequently, the least squares or weighted least squares method is used to solve for the long-interval three-dimensional relative position increment vector in the geocentric Earth-fixed coordinate system. Because the horizontal distance baseline between the reference epoch and the current epoch used in the long-interval difference method is much longer than that of adjacent epochs, the proportion of instantaneous turbulence displacement caused by sea waves in the total displacement of the long baseline is significantly diluted. Therefore, the long-interval relative position increment is mainly contributed by the ship's actual horizontal navigation displacement, exhibiting natural anti-turbulence characteristics. After obtaining the relative position increment in the geocentric-fixed coordinate system, the same coordinate transformation method as in the short-interval difference method is used to transform it to the ship's local northeast-sky coordinate system, obtaining the long-interval eastward displacement increment and the long-interval northward displacement increment. Finally, the heading of the above two horizontal displacement increments is calculated using the arctangent function, and quadrant correction is completed to obtain the long-interval coarse heading.

[0064] Through the above steps, the long-interval differential calculation process introduces a complete ephemeris orbital error compensation and atmospheric delay compensation stage before cross-epoch differential calculation. This effectively eliminates the systematic displacement deviation caused by broadcast ephemeris prediction errors and atmospheric delay changes under long timescale intervals. This allows for high-precision relative displacement calculation results even when the time interval between the reference epoch and the current epoch is significantly extended, thus ensuring the calculation accuracy of long-interval coarse headings and providing a highly stable heading to be verified for subsequent heading integrity testing.

[0065] Further, step S321 involves solving for the satellite orbit coordinates and satellite clock errors at two epochs, and calculating the line-of-sight orbit residual error, including: S3211 reads the satellite broadcast ephemeris corresponding to the current epoch and the reference epoch respectively, and calculates the satellite orbit coordinates and satellite clock difference at the two epochs.

[0066] For the current epoch and the long-interval reference epochs selected in step S3, the satellite broadcast ephemeris parameters received and decoded by the receiver at their respective epochs are obtained. Broadcast ephemeris is a set of orbit prediction data broadcast in real time by GNSS satellites. The receiver obtains Kepler orbital elements and their perturbation correction parameters, including the orbital semi-major axis, eccentricity, orbital inclination, right ascension of the ascending node, perigee distance, and mean perigee angle, as well as the polynomial fitting coefficients of the satellite clock offset, by demodulating the satellite navigation message. In the actual operation of a single-antenna GNSS receiver, broadcast ephemeris is usually updated at regular intervals, and the receiver internally stores the latest valid ephemeris data for each visible satellite. Based on the timestamps of the current epoch and the reference epochs, the corresponding valid broadcast ephemeris parameters are extracted from the receiver's ephemeris buffer. For possible ephemeris updates between two epochs, it must be ensured that each epoch uses the actual valid ephemeris version for its respective time. Using the extracted ephemeris parameters and the corresponding satellite orbit calculation models for each satellite navigation system, the three-dimensional orbital coordinates of each common-view satellite in the geocentric-geo-fixed coordinate system at two epochs, along with the corresponding satellite clock bias values, were calculated. Satellite clock bias refers to the deviation of the satellite's onboard atomic clock from the system's standard time, and is typically used in the form of equivalent distance for subsequent observation correction. The calculated satellite orbital coordinates and satellite clock biases constitute a set of spatiotemporally corresponding satellite state information, providing fundamental data for the subsequent quantitative calculation of orbital deviations and clock bias changes.

[0067] S3212 calculates the satellite orbital position deviation and satellite clock error difference between two epochs.

[0068] Using the satellite orbital coordinates and satellite clock difference obtained in step S3211 at the current epoch and reference epoch, the orbital position deviation vector and satellite clock difference difference are calculated for each common-view satellite with valid observations at both epochs. The orbital position deviation vector is calculated by directly taking the three-dimensional vector difference between the satellite orbital coordinates at the two epochs. The direction of this vector represents the direction of the satellite's displacement in space between the two epochs, and the magnitude of the vector represents the magnitude of the change in the satellite's orbital position. Because the orbital prediction accuracy of broadcast ephemeris is limited, typically on the order of meters to several meters, and the prediction error has the characteristic of accumulating over time, when the time interval between the current epoch and the reference epoch reaches several seconds or even longer, there will be a significant difference between the orbital position calculated from the broadcast ephemeris at the two epochs and the actual orbit of the same satellite. This difference is directly reflected in the inconsistency between the calculated orbital coordinates at the two epochs, i.e., the orbital position deviation vector. The satellite clock error difference is calculated by taking the algebraic difference between the satellite clock error values ​​at two epochs. Due to the limited frequency stability of the satellite's onboard atomic clock, the satellite clock error also changes slowly over time. The difference between the clock error values ​​at two moments represents the drift of the satellite clock error within that time interval. Orbital position deviation and satellite clock error difference together constitute the source of error introduced by broadcast ephemeris inconsistencies in long-interval differential calculations. If the original carrier phase observations are directly differentiald across epochs without compensation, the above errors will be fully propagated into the solution of the relative position increment, forming a systematic displacement deviation positively correlated with the time interval.

[0069] S3213 projects the satellite orbital position deviation and satellite clock error difference onto the satellite-receiver line-of-sight direction to obtain the line-of-sight orbital residual error.

[0070] The satellite orbital position deviation vector and satellite clock error difference calculated in step S3212 are projected along the line-of-sight direction to transform them into a one-dimensional equivalent distance error that directly affects the carrier phase observation. The carrier phase observation is essentially a distance measurement from the satellite to the receiver along the line-of-sight direction. The position deviation of the satellite orbit in three-dimensional space does not entirely affect this distance measurement; only the component along the line connecting the satellite and the receiver is perceived by the carrier phase observation. Therefore, it is necessary to project the three-dimensional orbital position deviation vector onto the satellite-receiver line-of-sight direction. This projection operation is achieved by calculating the dot product of the orbital position deviation vector and the line-of-sight direction unit vector. The line-of-sight direction unit vector is determined by the receiver's approximate position coordinates and the satellite's orbital coordinates at two epochs. For the current epoch and the reference epoch, their respective line-of-sight direction vectors can be calculated. During projection, the average of the line-of-sight directions at the two epochs can be taken, or the average can be taken after separate projections, to more accurately reflect the impact of the orbital deviation on the signal propagation path. The satellite clock error difference itself is a time quantity or equivalent distance quantity, and it is directly included in the compensation as the equivalent distance error along the line-of-sight direction. By combining the line-of-sight projection component of the orbital position deviation with the satellite clock error difference, the line-of-sight orbital residual error of each satellite between two epochs is obtained. This error term fully describes the systematic deviation caused by the inconsistency of broadcast ephemeris, which cannot be naturally eliminated in long-interval cross-epoch differences, and provides an accurate compensation basis for pre-correcting the carrier phase observations in subsequent steps.

[0071] Furthermore, in step S322, the ionospheric model and the tropospheric model are used to calculate the atmospheric delay difference between two epochs, including: S3221 uses the Klobuchar ionospheric model to calculate the ionospheric delay at the current epoch and the reference epoch, respectively, and obtains the ionospheric delay difference.

[0072] The Klobuchar ionospheric model was used to estimate the ionospheric delay of signal propagation at both the current epoch and a long-interval reference epoch. The Klobuchar model is an empirical ionospheric correction model widely used in GPS single-frequency receivers. This model simplifies the ionosphere to a thin shell approximately 350 km above the Earth's surface, assuming that all free electrons in the ionosphere are concentrated within this shell. The diurnal variation of the ionospheric vertical delay is described by a set of eight ionospheric parameters broadcast from the ephemeris. In the actual calculation, firstly, using the latitude and longitude coordinates obtained from the receiver's current epoch and the orbital coordinates of each visible satellite, the elevation and azimuth angles of the satellites relative to the receiver were calculated, thus determining the geographic coordinates of the puncture point of the signal propagation path on the ionospheric shell. Then, based on the geomagnetic latitude and local time of the puncture point, combined with the ionospheric alpha and beta parameters broadcast from the ephemeris, the ionospheric vertical delay along the propagation paths of each common-view satellite at both epochs was calculated. Finally, the vertical delay is mapped to an oblique ionospheric delay along the signal propagation direction using a tilt factor related to the satellite elevation angle. After independently calculating the ionospheric delay at two epochs, the ionospheric delay value at the current epoch is subtracted from the ionospheric delay value at the reference epoch to obtain the ionospheric delay difference for each common-view satellite between the two epochs. The variation in ionospheric delay mainly originates from the change in the distance the signal propagation path travels in the ionosphere caused by changes in satellite elevation angle, as well as the evolution of the total number of ionospheric electrons over time and space. In long-timescale cross-epoch difference calculations, this difference can reach the equivalent distance order of several centimeters to tens of decimeters. Without compensation, the accuracy of the relative displacement calculation will be directly degraded.

[0073] S3222 uses the Saastamoinen tropospheric model to calculate the tropospheric delay at the current epoch and the reference epoch, respectively, and obtains the tropospheric delay difference.

[0074] The Saastamoinen tropospheric model was used to estimate the tropospheric delay of the signal propagation at both the current epoch and a long-interval reference epoch. Tropospheric delay is the ranging error caused by the slowing of signal propagation speed and the curvature of the propagation path due to changes in atmospheric density as the GNSS signal passes through the Earth's neutral atmosphere. The Saastamoinen model is a classic model for calculating tropospheric zenith delay. This model divides the tropospheric delay into two components: a dry component and a wet component. The dry component, mainly caused by atmospheric hydrostatic equilibrium, accounts for about 90% of the total tropospheric delay and is relatively stable. The wet component, contributed by the water vapor content in the atmosphere, accounts for a smaller proportion but varies significantly in time and space. The input parameters required for the model calculation include meteorological elements such as surface atmospheric pressure, temperature, and water vapor partial pressure. In practical shipborne single-antenna GNSS applications, if the receiver is not equipped with a meteorological sensor, standard atmospheric parameters can be used as approximate inputs, such as a surface temperature of 15 degrees Celsius, an atmospheric pressure of 1013.25 hPa, and a typical value for the corresponding relative humidity for the water vapor partial pressure. After calculating the tropospheric zenith delay for each common-view satellite at two epochs, the zenith delay is projected onto the signal propagation direction using a mapping function related to the satellite elevation angle. Commonly used mapping functions include the Niell mapping function. After independently calculating the tropospheric delay at the two epochs, the tropospheric delay value at the current epoch is subtracted from the tropospheric delay value at the reference epoch to obtain the tropospheric delay difference for each common-view satellite between the two epochs. At long timescales, the tropospheric delay difference is also significant due to changes in satellite elevation angle and evolution of atmospheric temperature, humidity, and pressure conditions.

[0075] S3223 combines the ionospheric delay difference with the tropospheric delay difference to obtain the atmospheric delay difference.

[0076] The ionospheric delay difference calculated in step S3221 and the tropospheric delay difference calculated in step S3222 are combined to form the atmospheric delay difference used for pre-correction of carrier phase observations. The effects of ionospheric delay and tropospheric delay on GNSS signal propagation are physically independent. Ionospheric delay is mainly caused by the dispersion effect of free electrons on electromagnetic waves; its magnitude is inversely proportional to the square of the signal frequency, and its influence on carrier phase and pseudorange observations has opposite signs. Tropospheric delay, on the other hand, is independent of signal frequency and has the same sign on both carrier phase and pseudorange observations. Since the two types of atmospheric delays are not coupled in terms of physical causes and variation patterns, the ionospheric delay difference and tropospheric delay difference of the same satellite between two epochs can be directly algebraically added to obtain the corresponding atmospheric delay difference for that satellite. In practical processing, the atmospheric delay difference needs to be converted into a corresponding range correction based on the carrier phase frequency. For receivers using multi-frequency observations, the ionospheric delay at each frequency point can be scaled according to the frequency ratio. The merged atmospheric delay difference directly reflects the systematic bias of the carrier phase observations introduced by the difference in atmospheric conditions between two epochs in long-timescale cross-epoch difference calculations. This bias will be used together with the line-of-sight orbit residual error in subsequent step S323 to pre-correct the carrier phase observations, thereby eliminating the influence of atmospheric delay variations on the accuracy of relative position increment calculations before the difference operation.

[0077] Further, in step S324, cross-epoch difference is performed on the pre-corrected carrier phase observations to solve for the long-interval relative position increment, and the long-interval coarse heading is obtained through coordinate transformation and heading calculation, including: S3241 performs cross-epoch difference on the pre-corrected carrier phase observations to solve for the relative position increment in the long-interval geocentric-geofixed coordinate system.

[0078] Cross-epoch differential calculations are performed on the carrier phase observations after ephemeris orbital error compensation and atmospheric delay compensation to solve for the relative position increment between the current epoch and a long-interval reference epoch. After pre-correction, the long-timescale systematic error in the carrier phase observations has been significantly reduced, and the cross-epoch differential observation equations at this point mainly reflect the true spatial displacement relationship of the receiver between the two epochs. Specifically, the pre-corrected carrier phase observations of each common-view satellite at the current epoch and the reference epoch are subtracted between epochs to form differential observations. Since the pre-correction has compensated for the errors introduced by orbital deviation and atmospheric delay variations, the main parameter to be estimated in the differential observations is the three-dimensional position increment of the receiver in the geocentric-ground-fixed coordinate system between the two epochs. In the actual solution process, the observation weights can be set according to the carrier-to-noise ratio or satellite elevation angle of each satellite, assigning higher weights to satellites with better signal quality to suppress the influence of observation noise and multipath effects on the solution results. For satellites that may experience cycle slips between two epochs, cycle slip detection markers can be used to reduce their weight or remove them, ensuring the continuity of carrier phase observations used in the solution. Subsequently, least squares or weighted least squares parameter estimation methods are used to solve for the three-dimensional relative position increment vector in the geocentric-ground-fixed coordinate system from the differential observations of all common-view satellites. Since the reference epoch used in long-interval differential calculations has a large time span and spatial distance from the current epoch, the ship's actual navigation displacement within the long timescale dominates the total displacement, while the instantaneous turbulence displacement caused by sea waves is significantly diluted in the long baseline due to its zero-mean characteristic. Therefore, the long-interval relative position increment obtained has a natural anti-turbulence characteristic and can stably reflect the ship's macroscopic navigation direction.

[0079] S3242 transforms the relative position increment in the long-interval geocentric-ground-fixed coordinate system to the ship's local northeast-sky coordinate system, separates the vertical undulation displacement, and obtains the long-interval eastward displacement increment and the long-interval northward displacement increment.

[0080] The three-dimensional relative position increments obtained in the long-interval geocentric-fixed coordinate system obtained in step S3241 are transformed to the local northeast-sky coordinate system centered on the ship's current position. The coordinate transformation method is consistent with that used in step S22 of the short-interval difference process, that is, a rotation matrix from the geocentric-fixed coordinate system to the local northeast-sky coordinate system is constructed using the ship's latitude and longitude coordinates at the current epoch. Then, the long-interval three-dimensional relative position increment vector is multiplied by the rotation matrix to obtain three mutually orthogonal displacement components: eastward, northward, and celestial. Based on this, the celestial displacement component is separated and removed, retaining only the eastward and northward displacement components as long-interval horizontal displacement increments. In the long-interval difference process, the celestial displacement component also mainly comes from tidal level changes, average hull draft changes, and long-period heave motion caused by waves, and is unrelated to the ship's heading in the horizontal plane. After removing it, the long-interval eastward and northward displacement increments purely reflect the ship's cumulative displacement in the horizontal plane between the current epoch and the reference epoch, providing input quantities that are not disturbed by vertical motion for subsequent heading calculations. It is worth noting that, since the time baseline of long-interval difference is longer, the horizontal displacement increment is usually much larger than that of short-interval difference. Therefore, the significance of vertical separation is not only to eliminate interference, but also to ensure the geometric strength of the heading solution. The larger the horizontal displacement, the less sensitive the heading angle is to small disturbances, and the better the stability of the solution.

[0081] S3243 uses the arctangent function to calculate the heading of the long-interval eastward displacement increment and the long-interval northward displacement increment, thus obtaining the long-interval coarse heading.

[0082] The heading angle is calculated using the long-interval eastward and northward displacement increments obtained in step S3242. The basic method for calculating the heading angle is to use the northward displacement increment as the denominator and the eastward displacement increment as the numerator to calculate the arctangent value, obtaining the angle from the northward direction to the ship's horizontal direction of motion, rotating clockwise. Since the arctangent function typically ranges from -90 degrees to +90 degrees, covering only two quadrants, quadrant correction is required based on the signs of the eastward and northward displacement increments: when the northward displacement increment is positive, the heading angle is in the first or second quadrant; when the northward displacement increment is negative, the heading angle is in the third or fourth quadrant. Simultaneously, the sign of the eastward displacement increment is used to further distinguish between the left and right quadrants. After quadrant correction, the heading angle is normalized to the range of 0 to 360 degrees, with 0 degrees corresponding to true north, 90 degrees to true east, 180 degrees to true south, and 270 degrees to true west, consistent with the commonly used heading definition in navigation. This heading angle is the long-interval coarse heading. Since the long-interval coarse heading is calculated based on a large horizontal displacement, it has a significant suppression effect on high-frequency disturbances caused by instantaneous sea surface turbulence, and can stably indicate the average navigation direction of the ship over a longer time window. However, the reliability of the long-interval coarse heading depends on the quality of the carrier phase observations during the long-interval differential process and the adequacy of error compensation. In cases where satellite signals are blocked, multipath effects are severe, or the reference epoch is improperly selected, calculation anomalies may still occur. Therefore, it is necessary to perform subsequent heading integrity checks on it together with the short-interval reference heading.

[0083] Through the above steps, after error pre-correction, the long-interval differential calculation uses cross-epoch carrier phase differential to solve for the relative position increment with a long spatial baseline. After coordinate transformation and vertical separation, the pure horizontal displacement component is obtained. Finally, the long-interval coarse heading with high anti-turbulence capability is calculated, providing a heading value with significantly better stability than traditional single-epoch Doppler heading for subsequent heading integrity detection.

[0084] The inherent anti-turbulence properties of long-interval differential spectroscopy stem from the fundamental difference in statistical characteristics over time between the instantaneous turbulence displacement caused by sea waves and the ship's actual navigation displacement. The rolling and heaving motions caused by sea waves are high-frequency periodic disturbances; within a time window of several to tens of seconds, the direction of the resulting displacement vector changes continuously, with the statistical expectation approaching zero. In contrast, the actual horizontal displacement of a ship traveling at a constant speed has a definite direction and continuous accumulation, with the magnitude of the displacement vector increasing linearly with time. When using long-interval cross-epoch differential spectroscopy, the relative horizontal displacement between two epochs is composed of the ship's actual navigation displacement and the instantaneous turbulence displacement. As the time interval increases, the proportion of the actual navigation displacement in the total displacement continuously rises, while the contribution of the turbulence displacement decreases accordingly, significantly improving the signal-to-interference ratio. Taking a ship traveling at approximately 3 knots in sea state 4 as an example, the ship's actual horizontal displacement is about 15 meters within a 10-second time interval. While the instantaneous displacement amplitude caused by waves may reach several meters, its vector direction is variable, and its overall contribution to the long baseline is far less than the actual navigation displacement. Therefore, the heading angle calculated based on the relative displacement of the long baseline can stably reflect the ship's macroscopic navigation direction, effectively suppressing sea surface turbulence disturbances from a measurement mechanism perspective.

[0085] Specifically, in step S4, using the short-interval reference heading as the verification benchmark, a heading autonomy integrity test is performed on the long-interval coarse heading to obtain the stable heading output for the current epoch, including: S41, extract the short-interval reference heading from the short-interval buffer, obtain the long-interval coarse heading, and calculate the absolute value of the heading difference between the two.

[0086] The short-interval reference heading calculated in step S2 and stored in the short-interval buffer is numerically compared with the long-interval coarse heading calculated in step S3. The short-interval reference heading is a real-time heading obtained based on the phase difference of adjacent epoch carriers. Its time baseline is extremely short, and there is almost no error accumulation effect, which can reflect the horizontal motion direction perceived by the receiver at the current instant more realistically. The long-interval coarse heading is a stable heading obtained based on the difference of long baselines across epochs. Its advantage is that it dilutes the instantaneous disturbance caused by sea surface turbulence through long-distance relative displacement, but at the cost of being affected by abnormal factors such as cycle slips, multipath effects, or residual orbital errors over long time intervals. Although both are measurements of the same ship's heading at the current moment, their time scales and error characteristics are completely different. Therefore, the difference between the two can serve as an effective basis for judging whether the long-interval coarse heading is affected by abnormal interference. First, the short-interval reference heading value corresponding to the current epoch is read from the short-interval buffer, and the long-interval coarse heading value is obtained from the calculation result of step S3. Then, the absolute value of the difference between the two is calculated. Since the heading angle is an angular quantity that varies periodically within the range of 0 to 360 degrees, the difference between two heading angles needs to be calculated correctly based on the cyclical characteristics of the angle. For example, when two heading angles are close to 0 degrees and close to 360 degrees respectively, their actual angle difference should be taken as the value around the shorter arc, rather than simply being algebraically subtracted. The absolute value of the heading difference quantifies the degree of consistency between heading measurement results at two different time scales and is the core input for subsequent integrity assessment.

[0087] S42, compare the absolute value of the heading difference with a preset integrity threshold. If the absolute value of the heading difference is less than or equal to the preset integrity threshold, then the long-interval coarse heading is taken as the valid output heading.

[0088] The absolute value of the heading difference calculated in step S41 is judged using a preset heading integrity threshold. The preset integrity threshold is the boundary value for determining whether there are anomalies in the long-interval differencing results. Its selection needs to comprehensively consider the reasonable deviation range between the two heading solutions when the ship is sailing normally. When the ship is sailing in a straight line at a constant speed, both the short-interval reference heading and the long-interval coarse heading reflect the ship's direction of motion in the same horizontal plane. Since the short-interval heading contains high-frequency disturbance components while the long-interval heading is smoothed through a long baseline, a small deviation between the two is normal. This deviation is usually within the range allowed by the preset integrity threshold. When the absolute value of the heading difference is less than or equal to the preset integrity threshold, it indicates that the long-interval differencing results have good consistency with the short-interval reference heading, the long-interval differencing process is not significantly interfered with by abnormal factors such as cycle slips, severe multipath effects, or residual track errors, and the solution results are reliable. At this point, the long-interval coarse heading is determined as the effective output heading that has passed the integrity test. Its advantage lies in the fact that the long baseline displacement has a natural suppressive effect on sea surface turbulence, and the output heading has higher stability, which can provide smooth and continuous heading reference information for the ship control system.

[0089] S43, if the absolute value of the heading difference is greater than the preset integrity threshold, it is determined that there is an anomaly in the long-interval coarse heading, and the short-interval reference heading is output as the alternative heading for the current epoch.

[0090] When the absolute value of the heading difference exceeds a threshold, it indicates an out-of-normal directional deviation between the long-interval differential result and the short-interval reference heading, casting doubt on the reliability of the long-interval coarse heading. Possible causes of this anomaly include, but are not limited to, the following: cycle slips occur in the carrier phase of some common-view satellites in the long-interval reference epoch or the current epoch, introducing integer cycle slip errors into the differential observation equations; during the long-interval differential process, despite compensation for orbital errors and atmospheric delays, the residual system bias still exceeds the expected range; the observations of one or more satellites are severely affected by multipath effects, such as signal reflections caused by nearby dock structures or large vessels when a ship is berthing; or, due to limited buffer time, a reference epoch with insufficient spatial distance is selected during reference epoch selection, resulting in insufficient signal-to-noise ratio of the long baseline displacement itself. In these situations, continuing to use the long-interval coarse heading as the output may lead to erroneous heading information being transmitted to the ship's control system, affecting navigation safety. Therefore, the heading integrity warning mechanism is triggered, abandoning the long-interval coarse heading of the current epoch and outputting the short-interval reference heading as an alternative heading. Although the short-interval reference heading contains a high instantaneous disturbance component under turbulent sea conditions, it is calculated based on an extremely short time baseline, eliminating the risk of error accumulation and anomaly propagation. It can reliably reflect the current instantaneous direction of motion and provides effective redundancy protection when anomalies occur in the long-interval differential heading.

[0091] S44. Perform quadrant correction and normalization on the output heading determined in step S42 or step S43 to obtain stable heading data for the current epoch.

[0092] Regardless of whether the heading output after the integrity check comes from the long-interval coarse heading in step S42 or the short-interval reference heading in step S43, it undergoes final quadrant correction and normalization to ensure the uniformity and standardization of the output heading data format. Quadrant correction aims to ensure that the heading angle is correctly mapped to a 0-360 degree coordinate system with true north as the 0-degree reference, eliminating quadrant ambiguities that may arise due to the limitations of the arctangent function's range. Specifically, based on the signs of the eastward and northward displacement components constituting the heading, the original calculated angle values ​​are corrected to the correct quadrant. Normalization aims to standardize the numerical range of the heading angle. If the angle value after quadrant correction is still outside the 0-360 degree range, or if negative values ​​or values ​​greater than 360 degrees appear in certain boundary cases, it is adjusted to the standardized range by adding or subtracting an integer period of 360 degrees. After quadrant correction and normalization, stable heading data for the current epoch is obtained. This data is output in the common true north reference heading format used in the navigation field and can be directly read and used by the ship's heading display instruments, autopilot, or navigation management system.

[0093] Through the above four steps, the autonomous heading integrity detection process uses the short-interval reference heading as the redundancy verification benchmark, performs consistency comparison and anomaly detection on the long-interval coarse heading, prioritizes the output of its high-stability heading when the long-interval differential result is normal, and automatically switches to the short-interval reference heading as a substitute when an anomaly is detected. Quadrant correction and normalization are uniformly performed before output, thus achieving autonomous heading integrity monitoring and seamless switching under fault conditions without the need for any external reference information, while taking into account the high stability and high reliability of heading output.

[0094] Sea trials have verified that, under sea state 4 and a low-speed, uniform navigation condition of 2 to 5 knots, the traditional heading calculation scheme based on single-epoch Doppler instantaneous velocity vectors exhibits heading jitter of ±8 to ±15 degrees. However, the heading stability measurement method provided in this embodiment effectively controls heading jitter within ±1.2 degrees, improving heading stability by over 85% compared to the traditional method, with no significant output lag. It should be noted that the above test data are exemplary results under specific sea states and navigation conditions and do not constitute a limitation on the scope of protection of this invention.

[0095] Furthermore, in one embodiment of the present invention, after obtaining the stable heading output for the current epoch in step S4, the method further includes: S51, obtain the speed and heading values ​​for a consecutive preset number of observation epochs, and calculate the speed fluctuation value and heading change rate.

[0096] The ship's navigation status is monitored in real time using the stable heading data output from multiple consecutive observation epochs in step S4, along with the corresponding speed information. The heading value is directly taken from the stable heading data output from each observation epoch in step S4. The speed value can be calculated using the ratio of the horizontal displacement increment obtained by short-interval differential calculation in step S2 to the corresponding epoch time interval. Specifically, the horizontal displacement increments of the short-interval eastward and northward directions are combined and then divided by the time difference between adjacent epochs to obtain the instantaneous speed estimate for each observation epoch. Alternatively, the long-timescale average speed can be obtained by using the ratio of the horizontal displacement increment obtained by long-interval differential calculation in step S3 to the time interval between two epochs. Under uniform straight-line navigation conditions, the two speed estimates show good consistency. A sliding time window is formed by a preset number of observation epochs. The window length can be set according to the receiver's sampling frequency and the ship's maneuvering response time requirements. For example, with a 10Hz sampling frequency, 5 to 10 consecutive epochs can correspond to an observation duration of 0.5 to 1 second. This captures the trend of motion changes within a short period without introducing significant decision delays due to an excessively long window. After acquiring the speed and heading values ​​at each epoch within the sliding window, two indicators characterizing the ship's motion stability are calculated: speed fluctuation and heading change rate. Speed ​​fluctuation can be measured using the standard deviation or maximum deviation of the speed sequence within the window, reflecting the dispersion of the ship's speed over a short period. The heading change rate can be measured using the average of the heading differences between adjacent epochs within the window or the ratio of the heading difference between the first and last epochs to the time interval, reflecting the rate of change of the ship's direction over a short period. These two indicators characterize the ship's current motion state from the dimensions of speed magnitude and angular velocity, respectively, providing a quantitative basis for subsequent navigation condition determination.

[0097] S52, based on the comparison results of the speed fluctuation value and the rate of change of course with the preset speed fluctuation threshold and the preset rate of change of course threshold, determine the current navigation condition of the ship.

[0098] The speed fluctuation value and heading change rate calculated in step S51 are compared with preset thresholds to determine the ship's current navigation condition. The preset speed fluctuation threshold and preset heading change rate threshold are empirical parameters set based on the ship's typical motion characteristics and navigational experience. For example, the speed fluctuation threshold can be set to 1.0 knots, and the heading change rate threshold can be set to 0.5 degrees per second. In the actual determination process, the speed fluctuation value is compared with the preset speed fluctuation threshold, and the heading change rate is compared with the preset heading change rate threshold. When the speed fluctuation value is less than the preset speed fluctuation threshold and the heading change rate is less than the preset heading change rate threshold, it indicates that the ship's speed and direction remain stable within the current slip time window, conforming to the motion characteristics of uniform straight-line navigation. In this case, the ship's current navigation condition is determined to be uniform straight-line navigation. This condition is common during the constant-speed cruising phase on ocean routes or the low-speed, stable berthing approach phase within port areas. Conversely, when the speed fluctuation value is greater than or equal to the preset speed fluctuation threshold, or the heading change rate is greater than or equal to the preset heading change rate threshold, it indicates that the ship's speed or direction has changed significantly, and the ship is performing acceleration, deceleration, turning, or combined maneuvering. In this case, the ship's current navigation condition is determined to be a turning or speed change maneuvering condition. This dual-indicator parallel judgment method for condition identification has low computational load and fast response speed. It can complete the classification and determination of navigation status in real time using only GNSS measurement data without relying on external sensors and electronic chart information.

[0099] In this embodiment, the determination of the uniform straight-line navigation condition in step S52 adopts a sliding window confirmation mechanism based on three consecutive observation epochs. Specifically, the current navigation condition is determined to be uniform straight-line navigation only when the speed fluctuation values ​​of three consecutive observation epochs are all less than a preset speed fluctuation threshold, and the heading change rate of three consecutive observation epochs is all less than a preset heading change rate threshold. This continuous confirmation mechanism can effectively avoid misjudgment of the navigation condition due to instantaneous measurement noise of a single epoch, and prevent output heading oscillation caused by frequent switching between steady-state and maneuver modes of the differential strategy.

[0100] S53 adaptively adjusts the differential strategy and detection threshold based on the determined navigation conditions.

[0101] Based on the navigation condition category determined in step S52, the solution strategy for cross-epoch difference and the relevant thresholds for heading integrity detection are dynamically adjusted to achieve an optimal balance between stability and response speed in the heading measurement method. Under uniform straight-line navigation conditions, the ship's direction of motion remains stable over a relatively long period, making the stability of the heading output the primary objective. Therefore, maintaining the complete long-short double-layer difference and heading integrity detection process from steps S2 to S4, long-interval coarse headings are prioritized as the main output heading, fully leveraging the suppression effect of long-baseline difference on sea surface turbulence disturbances to provide smooth and stable heading information for the ship's control system. Under turning or variable-speed maneuvering conditions, the ship's direction of motion is changing rapidly, making real-time heading output the primary objective. However, due to the large reference epoch time span, the output heading of long-interval difference naturally lags behind the ship's current true heading in the time domain. Continuing to prioritize long-interval results would lead to sluggish control response. Therefore, under this operating condition, the differential strategy is adaptively adjusted, including raising the preset time threshold and preset horizontal distance threshold to shorten the allowable time difference upper limit for long-interval reference epochs, raising the preset integrity threshold to relax the judgment criteria for the consistency of long-interval results, and reducing the output weight of long-interval coarse heading while increasing the output weight of short-interval reference heading. After the above adjustments, the contribution ratio of short-interval differential to heading output increases, the heading response speed is improved, and it can track the ship's maneuvers in a timely manner. After the ship returns to uniform straight-line navigation, the strategy automatically reverts to a steady-state mode dominated by long-interval differential. This mechanism of real-time adaptive switching based on operating conditions enables the heading measurement method to provide heading output performance that matches the navigation state under all operating conditions.

[0102] Furthermore, in step S53, based on the determined navigation conditions, the differential strategy and detection threshold are adaptively adjusted, including: S531, if it is determined that the ship is in a uniform straight-line navigation condition, continue to execute steps S2 to S4 in a loop, and use the long-interval coarse heading as the main output heading.

[0103] When step S52 determines that the ship is currently in a uniform speed straight-line navigation condition, the heading measurement method maintains the standard two-layer differential and integrity detection process unchanged. Uniform speed straight-line navigation is a common stable motion state for ships during ocean cruising, straight-line passages in waterways, or low-speed berthing in ports. In this state, the ship's heading remains basically constant over a relatively long time window, without significant changes in direction in a short period. At this time, the core requirement for heading output is to minimize heading jitter caused by sea surface turbulence, providing smooth and continuous heading reference information for the ship's control system and crew, enabling accurate track maintenance and heading monitoring. Based on this requirement, a strategy of using long-interval coarse headings as the primary output heading is adopted. Long-interval coarse headings are calculated by solving the relative displacement of long baselines between the ship and historical reference epochs. Due to the long baseline distance, the proportion of the ship's instantaneous rolling and heaving motions under wave action in the total displacement vector is significantly diluted, and the calculation results can stably reflect the ship's macroscopic navigation direction, with the output heading sequence exhibiting good smoothness and continuity. Meanwhile, the complete process from steps S2 to S4 continues to be executed in a loop. The short-interval reference heading is still calculated at each epoch and used as a verification benchmark to perform autonomous heading integrity checks on the long-interval coarse heading. This ensures that if an anomaly occurs in the long-interval differential during uniform speed navigation, the heading output can be switched to the short-interval heading output in a timely manner through the integrity check mechanism. This ensures that the reliability of the heading output is not sacrificed while pursuing high stability.

[0104] S532 If it is determined that the ship has entered a turning or speed change maneuver, the preset time threshold, preset horizontal distance threshold and preset integrity threshold are increased, and the maximum time difference of the reference epochs that can be screened in the long interval buffer is shortened.

[0105] When step S52 determines that the ship has entered a turning or speed-changing maneuver, the threshold parameters involved in the long-interval differential are adaptively increased. During a turning maneuver, the ship's heading angle continuously changes, potentially traversing tens of degrees or even larger angles within a short period. During a speed change, the ship's speed amplitude changes significantly, and its track curvature and motion state differ markedly from those during the uniform straight-line phase. Under such maneuvering conditions, if the historical reference epochs used in the long-interval differential have a large time span, the corresponding ship position and heading will differ significantly from the actual state at the current moment. The calculated long-interval coarse heading will then fail to accurately reflect the ship's current true heading, resulting in significant output lag. To address this issue, the preset time threshold and preset horizontal distance threshold are increased, making it easier to match historical epochs with shorter time intervals to the current epoch when filtering reference epochs in the long-interval buffer, thereby shortening the time difference between the actual reference epoch and the current epoch. Meanwhile, the maximum time difference of reference epochs allowed for filtering in the long-interval buffer is shortened, limiting the selection range of reference epochs at the buffer depth level and further reducing the possibility of outdated historical data being selected. Furthermore, under maneuvering conditions, since the ship's course is actually changing, the difference between the short-interval reference course and the long-interval coarse course includes not only potential abnormal deviations but also course tracking differences due to different time scales, thus expanding the normal deviation range compared to the uniform speed condition. Therefore, the preset integrity threshold is simultaneously increased, and the consistency judgment standard for the long-interval coarse course is appropriately relaxed to avoid reasonable course differences caused by normal ship turning or speed changes being misjudged as long-interval differential anomalies, ensuring the correct operation of the integrity detection mechanism under maneuvering conditions.

[0106] In step S532, the threshold can be adjusted proportionally or in fixed steps. When adjusting proportionally, the preset time threshold, preset horizontal distance threshold, and preset integrity threshold can be increased to 1.5 to 2 times the base value to accommodate reasonable deviations caused by rapid changes in heading during maneuvers. When adjusting in fixed steps, the integrity threshold can be increased by 2 to 3 degrees each time. Shortening the maximum time difference of the reference epochs that can be filtered in the long-interval buffer can be achieved by limiting the traversal to only the most recent proportion of historical data in the buffer, for example, searching only historical epochs within the last 10 seconds.

[0107] S533 reduces the output weight of long-interval coarse heading and increases the output weight of short-interval reference heading during steering or gear shifting maneuvers.

[0108] To address the real-time requirements of heading output during steering or gear-changing maneuvers, the contribution ratios of long-interval coarse headings and short-interval reference headings in the final output heading are dynamically adjusted. The inherent advantage of long-interval coarse headings lies in their resistance to turbulence and high stability; however, their calculation based on a long historical baseline means they essentially reflect the ship's average heading over a past period. In maneuvering scenarios with rapidly changing headings, this averaging effect becomes a disadvantage, causing the output heading to lag behind the actual heading. Short-interval reference headings, based on differences between adjacent epochs, have an extremely short time baseline. Although they contain more high-frequency disturbance components, they can respond to instantaneous changes in the ship's heading with almost no delay. During maneuvers, heading response speed takes precedence over heading smoothness; the helmsman needs to perceive the ship's ongoing directional changes promptly, rather than a lagging heading that has been smoothed over a long period. Therefore, by reducing the output weight of long-interval coarse headings and increasing the output weight of short-interval reference headings, the final output heading value is more significantly contributed by the short-interval difference results. The weighting can be adjusted gradually, with the weight coefficients switching progressively over several epochs after the maneuver condition determination takes effect, to avoid sudden changes in output heading caused by abrupt weight changes. After this adjustment, the dynamic response characteristics of the heading output are significantly improved, enabling close tracking of real-time heading changes during steering or speed changes, providing accurate and timely control feedback to the helmsman. Once the ship completes the maneuver and returns to uniform straight-line navigation, the strategy in step S531 re-activates, and the weight coefficients gradually revert to a steady-state mode dominated by long-interval, coarse headings, thus achieving an optimal balance between heading stability and response speed under all operating conditions.

[0109] In step S533, the adjustment of the output weights adopts a gradual transition method, gradually switching the weight coefficients over several epochs after the maneuver condition determination takes effect. For example, in the first epoch after the determination of entering the maneuver condition, the weight of the long-interval coarse heading can be reduced from 1.0 to 0.7, in the second epoch to 0.3, and in the third epoch to 0.0. Correspondingly, the weight of the short-interval reference heading can be gradually increased from 0.0 to 1.0 to avoid output heading jumps caused by sudden weight changes. When the ship returns to the uniform straight-line navigation condition, the weight coefficients gradually revert to the steady-state configuration dominated by the long-interval coarse heading in the same gradual manner.

[0110] Compared to existing heading calculation schemes based on Doppler instantaneous velocity vectors, this invention no longer relies on velocity measurements from a single epoch. Instead, it characterizes the ship's macroscopic navigation direction by solving for the relative displacement vector with a long spatial baseline through long-interval, cross-epoch differential calculations. This avoids the amplification effect of instantaneous velocity disturbances on heading calculations from a measurement mechanism perspective. Compared to existing improved schemes based on filtering smoothing or inertial sensor fusion, this invention avoids the output lag problem introduced by filtering and does not rely on additional hardware such as IMUs or electronic compasses. It achieves a pure software-level improvement in heading stability while ensuring real-time heading response. More importantly, existing single-antenna direction finding schemes do not provide heading integrity monitoring methods. This invention, through an autonomous heading integrity detection mechanism using short-interval reference heading as a verification benchmark, achieves autonomous identification and reliable switching of heading anomalies for the first time in a single-antenna architecture, filling a technological gap in this field.

[0111] The following provides three complete embodiments, corresponding to three typical application scenarios: low-speed berthing, constant-speed cruising, and maneuvering. Each embodiment is described according to the complete process of steps S1 to S5.

[0112] Example 1: Low-speed berthing of a ship This embodiment describes a course measurement scenario of a small to medium-sized cargo ship during low-speed berthing in a port. The ship is equipped with a single-antenna GNSS receiver with a sampling frequency of 10Hz. At the time of berthing, the ship's speed is approximately 2 to 3 knots, the sea state is level 4, and there is significant wind and wave turbulence.

[0113] In step S1, the receiver continuously acquires raw observation data at a frequency of 10Hz, including timestamps, pseudorange observations, carrier phase observations, carrier-to-noise ratio, and latitude and longitude coordinates for each epoch. Two independent sliding buffer queues are constructed. The long-interval buffer has a buffer duration of 30 seconds and can store approximately 300 historical epoch data points. The short-interval buffer retains only the data for the current epoch and the previous adjacent epoch, and is updated using a first-in, first-out (FIFO) method. Since the ship is in a low-speed berthing state, the horizontal displacement accumulates slowly, and the historical data over a longer time range in the long-interval buffer is crucial for selecting reference epochs that meet the spatial distance conditions.

[0114] In step S2, two sets of carrier phase observations—one for the current epoch and one for the previous adjacent epoch—are extracted from the short-interval buffer. First, the common-view satellites of the two epochs are matched, eliminating satellites with observations only in a single epoch. For the commonly observed satellites, epoch-level carrier phase differential observations are constructed separately. Since the interval between adjacent epochs is only 0.1 seconds, the carrier phase integer ambiguity remains constant under continuous signal tracking conditions and is automatically eliminated by subtracting between epochs. The least squares method is used to solve for the three-dimensional relative position increment vector in the geocentric-fixed coordinate system. Subsequently, a rotation matrix is ​​constructed using the ship's latitude and longitude coordinates of the current epoch, transforming the relative position increment to the ship's local northeast-sky coordinate system, obtaining three displacement components: eastward, northward, and celestial. During berthing, the ship experiences significant heave motion due to wave influence, with the celestial displacement component potentially reaching tens of decimeters; this component is discarded, retaining only the eastward and northward horizontal displacement components. After solving using the arctangent function and completing quadrant correction, the short-interval reference heading is obtained; for example, the angle calculated in the current epoch is 35.2 degrees. This heading is stored in the short-interval buffer as a verification reference.

[0115] In step S3, reference epochs are selected from the long-interval buffer. Historical epochs in the buffer are traversed, and the horizontal relative distance and time difference between each historical epoch and the current epoch are calculated. Since the ship's speed is only about 2 knots, the horizontal displacement is about 0.1 meters in 0.1 seconds. Reaching the preset 10-meter horizontal distance threshold would require approximately 100 seconds, far exceeding the 30-second buffer duration. Therefore, none of the historical epochs in the buffer meet the horizontal distance threshold condition. However, the time difference between the farthest historical epoch in the buffer and the current epoch reaches 30 seconds, satisfying the preset 5-second time threshold condition. The historical epoch with the largest time difference is selected as the reference epoch. Subsequently, the satellite broadcast ephemeris at the current epoch and the reference epoch are read, and the orbital coordinates and satellite clock differences of the common-view satellites at the two times are calculated. The orbital position deviation vector and satellite clock error difference between two epochs are calculated. The orbital deviation is projected onto the line-of-sight direction between the satellite and the receiver to obtain the line-of-sight orbital residual error. Simultaneously, the Klobuchar and Saastamoinen models are used to calculate the ionospheric delay difference and tropospheric delay difference. After pre-correcting the carrier phase observations using the above error terms, cross-epoch finite difference is performed to solve for the long-interval relative position increment in the geocentric Earth-fixed coordinate system. After coordinate transformation and vertical separation, the long-interval eastward and northward displacement increments are obtained, and the long-interval coarse heading is calculated, for example, 35.8 degrees.

[0116] In step S4, the short-interval reference heading of 35.2 degrees and the long-interval coarse heading of 35.8 degrees are extracted, and the absolute value of the difference between the two headings is calculated to be 0.6 degrees. With a preset integrity threshold of 2 degrees, 0.6 degrees is less than 2 degrees, so the long-interval coarse heading is deemed to have passed the integrity check and is taken as a valid output heading. After quadrant correction and normalization, the stable heading data for the current epoch, 35.8 degrees, is output. Because the long-interval differential is based on the calculation of approximately 30 meters of actual navigation displacement accumulated over 30 seconds, the instantaneous turbulence displacement caused by sea waves is significantly diluted in the total displacement. The jitter amplitude of the output heading can be controlled within ±1.2 degrees, which is significantly better than the jitter level of ±8 to ±15 degrees under the same operating conditions of traditional Doppler velocity measurement schemes.

[0117] In step S5, the speed and heading values ​​for three consecutive epochs are obtained. The speed values ​​are obtained by removing the short-interval horizontal alignment by the epoch time interval. The speeds for the three epochs are 2.1 knots, 1.9 knots, and 2.0 knots, with a speed fluctuation of approximately 0.1 knots. The heading values ​​are 35.8 degrees, 35.6 degrees, and 35.9 degrees, with a heading change rate of approximately 0.15 degrees per second. The speed fluctuation is less than a preset threshold of 1.0 knots, and the heading change rate is less than a preset threshold of 0.5 degrees per second. It is determined that the ship is in a uniform straight-line navigation condition, and the standard procedures of steps S2 to S4 are continuously executed in a loop.

[0118] Example 2: Ship cruising at a constant speed This embodiment describes a course measurement scenario of a fishing vessel cruising at a constant speed in a straight line in open water. The vessel is equipped with a single-antenna GNSS receiver with a sampling frequency of 10Hz, a cruising speed of approximately 8 knots, a sea state of 3, and moderate swells on the sea surface.

[0119] In step S1, the receiver continuously acquires raw observation data at a frequency of 10Hz, constructing a long-interval buffer and a short-interval buffer. The buffer duration of the long-interval buffer is set to 30 seconds. Due to the high speed of the ship, the horizontal displacement accumulates rapidly, providing favorable conditions for subsequent selection of historical epochs that meet the spatial distance requirements.

[0120] In step S2, differential calculations are performed on the carrier phase observations of adjacent epochs in the short-interval buffer to solve for the relative position increment in the geocentric Earth-fixed coordinate system. After coordinate transformation to the Northeast-Sky coordinate system, the celestial displacement component is removed to obtain the short-interval eastward and northward displacement increments. At a speed of 8 knots, the horizontal displacement within 0.1 seconds is approximately 0.41 meters. The short-interval reference heading, for example, 120.3 degrees, is obtained through arctangent calculation and quadrant correction. This reference heading is then stored in the short-interval buffer.

[0121] In step S3, the long-interval buffer is traversed to filter reference epochs. At a speed of 8 knots, the ship can accumulate approximately 10 meters of horizontal displacement within 2.5 seconds, therefore, there are multiple historical epochs in the buffer that meet the horizontal distance threshold condition. Using the current epoch as a reference, candidate epochs with a horizontal relative distance greater than or equal to 10 meters or a time difference greater than or equal to 5 seconds are selected. Among the multiple candidate epochs, the epoch with the shortest time interval from the current epoch is selected as the reference epoch, for example, a historical epoch about 2.5 seconds ago, at which point the horizontal distance has just reached the 10-meter threshold, and the time difference is approximately 2.5 seconds. Due to the short time interval, the accumulated amount of broadcast ephemeris orbital error and atmospheric delay is relatively small, which is beneficial to improving the accuracy of subsequent compensation. After completing the ephemeris orbital error compensation and atmospheric delay compensation, cross-epoch carrier phase differential is performed to solve for the long-interval relative position increment. Within 2.5 seconds, the ship's horizontal displacement is approximately 10.3 meters, and the long-interval coarse heading is calculated, for example, 120.1 degrees.

[0122] In step S4, the short-interval reference heading of 120.3 degrees and the long-interval coarse heading of 120.1 degrees are extracted. The absolute value of the heading difference is 0.2 degrees, which is much smaller than the preset integrity threshold of 2 degrees. The long-interval coarse heading is determined to be valid, and 120.1 degrees is output as the stable heading for the current epoch. Since the baseline length used in the long-interval differential is about 10 meters, the instantaneous turbulence displacement caused by the surge accounts for less than a few percent of the total displacement, resulting in good smoothness and continuity in the output heading.

[0123] In step S5, the speed and heading data are continuously monitored for three epochs. Since both speed fluctuations and heading change rates are below preset thresholds, the system is continuously classified as uniform straight-line navigation, maintaining a steady-state mode with long-interval coarse headings as the primary output. Throughout the cruise, the heading output fluctuates stably around 120 degrees without significant jitter.

[0124] Example 3: Ship Maneuvering and Turning Condition This embodiment describes a course measurement scenario of a passenger ship performing a turning maneuver in a port channel. The ship is equipped with a single-antenna GNSS receiver with a sampling frequency of 10Hz, an initial speed of approximately 5 knots, and a course that gradually changes from 45 degrees to 135 degrees during the turning process, with a course change rate of approximately 3 degrees per second.

[0125] In steps S1 and S2, the receiver normally acquires data and constructs a two-layer buffer, calculating the short-interval reference heading epoch-by-epoch using short-interval differential calculation. Before the turn begins, the ship is in a constant-speed, straight-line state, and the short-interval reference heading is stable at around 45 degrees. After the turn begins, the short-interval reference heading is updated epoch-by-epoch with the ship's turning maneuvers: 45.3 degrees at the first epoch, 48.1 degrees at the second epoch, and 51.2 degrees at the third epoch, tracking the ship's heading changes in real time.

[0126] In step S3, the long-interval reference epoch selection and differential process continues. Since the ship maintained a constant speed and straight course for approximately 20 seconds before turning, the long-interval buffer contains abundant historical data. The horizontal distance between the current epoch and the historical epoch approximately 2.5 seconds ago is approximately 6.4 meters, which has not yet reached the 10-meter horizontal distance threshold, but the time difference has not yet reached 5 seconds. At this point, the historical epoch approximately 5 seconds ago might be selected as the reference epoch, at which point the horizontal distance is approximately 12.8 meters, satisfying the horizontal distance threshold condition. However, the ship's heading 5 seconds ago was still approximately 45 degrees, while the ship's heading at the current epoch has changed to approximately 51 degrees. The long-interval coarse heading calculation reflects the average direction between the two epochs, which may be approximately 48 degrees, showing a significant lag compared to the current actual heading.

[0127] In step S4, the short-interval reference heading is 51.2 degrees, and the long-interval coarse heading is 48.0 degrees, with an absolute difference of 3.2 degrees. Since the ship's heading is actually changing, this 3.2-degree difference includes both potential abnormal deviations and differences in heading tracking due to different time scales. If the integrity threshold remains at 2 degrees, this difference will trigger an anomaly detection. However, at this point, the operational condition identification mechanism in step S5 has already intervened.

[0128] In step S5, speed and heading data for three consecutive epochs are acquired. The heading values ​​for the three epochs are 45.3 degrees, 48.1 degrees, and 51.2 degrees, respectively, with a heading change rate of approximately 3 degrees per second, far exceeding the preset threshold of 0.5 degrees per second, indicating that the ship has entered a maneuvering condition. Based on the adaptive adjustment strategy under the maneuvering condition, the preset integrity threshold is increased from 2 degrees to 5 degrees, and the preset time threshold and horizontal distance threshold are simultaneously increased. The allowable time difference limit for reference epochs in the long-interval buffer is shortened to within the most recent 10 seconds. Simultaneously, the output weight of the long-interval coarse heading is reduced, while the output weight of the short-interval reference heading is increased. After the threshold adjustment, the heading difference of 3.2 degrees is less than the new integrity threshold of 5 degrees, and the long-interval coarse heading is not considered abnormal. However, due to the weight adjustment, the final output heading reflects more the real-time changes of the short-interval reference heading; for example, the overall output is approximately 50.5 degrees. Once the ship completes the turn and resumes a uniform straight course of 135 degrees, the rate of change of course drops below the threshold for three consecutive epochs, indicating a return to uniform speed. The integrity threshold and other parameters automatically revert to their base values, and the output weights are tilted back towards the long-interval coarse course, gradually making the course output smoother and more stable.

[0129] As can be seen from the above three embodiments, this measurement method can automatically adapt the optimal differential strategy and detection parameters according to the actual navigation state of the ship. In the scenario of strong turbulence during low-speed berthing, it can give full play to the anti-disturbance advantage of long-interval differential, continuously output high-stability heading in the scenario of uniform speed cruise, and improve the heading response speed in the scenario of maneuvering and turning. It achieves the optimization of heading measurement performance under all working conditions.

[0130] Accordingly, a second aspect of the present invention provides an electronic device, comprising: at least one processor; and a memory connected to the at least one processor; wherein the memory stores instructions executable by the at least one processor, the instructions being executed by the at least one processor to cause the at least one processor to perform the above-described method for measuring the anti-turbulence heading stability of marine GNSS based on cross-epoch difference.

[0131] Accordingly, a third aspect of the present invention provides a computer-readable storage medium having computer instructions stored thereon, which, when executed by a processor, implement the above-described method for measuring the anti-turbulence heading stability of marine GNSS based on cross-epoch difference.

[0132] The embodiments of this invention aim to protect a method for measuring the anti-turbulence heading stability of marine GNSS based on cross-epoch difference, which has the following effects: 1. By constructing a two-layer architecture that coordinates short-interval differential and long-interval cross-epoch differential operations, and using the short-interval differential heading as the verification benchmark, this system innovatively achieves receiver autonomous integrity monitoring in the heading dimension. The short-interval differential automatically eliminates integer ambiguity by utilizing the carrier phase difference between adjacent epochs, and can output a reference heading without accumulated error in real time without relying on external sensors. The long-interval differential, by selecting historical epochs that meet time or distance thresholds, uses the long baseline relative displacement vector to characterize the ship's macroscopic navigation direction, which significantly dilutes the instantaneous position and velocity disturbance errors introduced by sea surface wind and wave turbulence from a mechanistic perspective. On this basis, by comparing the consistency between the short-interval reference heading and the long-interval coarse heading, when the long-interval result is abnormal due to cycle slip, multipath, or residual orbit error, it can autonomously determine and switch to the short-interval heading output, fundamentally solving the defect of traditional single-antenna direction finding schemes lacking integrity monitoring capabilities at the heading level, and significantly improving the reliability of heading output. 2. To address the technical challenges of systematic displacement deviations caused by the accumulation of satellite broadcast ephemeris orbital bias, satellite clock bias, and ionospheric and tropospheric delays over long time in long-time-difference cross-epoch differential scenarios, this paper proposes a method that proactively reads the broadcast ephemeris corresponding to the current epoch and the reference epoch before performing long-interval cross-epoch differential calculations. This method calculates the satellite orbital coordinates and satellite clock bias at both times, calculates the line-of-sight orbital residual error, and simultaneously compensates for the atmospheric delay difference between the two epochs using ionospheric and tropospheric models. After pre-correcting the carrier phase observations using all the aforementioned error terms, the differential calculation is performed. This compensation process effectively eliminates the impact of long-timescale systematic displacement shifts on relative position increments, improves the upper limit of the calculation accuracy of long-interval differentials, and enables high-precision heading measurement results even with significantly extended time intervals, thus broadening the applicable time range for cross-epoch differential calculations. 3. By monitoring the ship's navigation status in real time, the system automatically identifies uniform straight-line navigation and turning / variable-speed maneuvering conditions, and adaptively adjusts the differential strategy and detection thresholds according to the condition. During uniform cruising or low-speed berthing, long-interval differential results are prioritized as the primary output heading, fully leveraging their advantages in anti-turbulence and high stability. During turning or variable-speed maneuvering, the system automatically increases the time threshold, distance threshold, and integrity detection threshold, shortens the allowable time difference upper limit of long-interval reference epochs, reduces the weight of long-interval differential output, and increases the weight of short-interval differential output, thus ensuring both heading response speed and output stability. This adaptive mechanism resolves the inherent contradiction that fixed differential strategies cannot simultaneously satisfy high heading stability and high real-time performance, enabling this solution to adapt to the heading measurement needs of ships in various sea states and all operating conditions, possessing universal capability across all scenarios.

[0133] 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 embodied 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.

[0134] This application is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of this application. 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... Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0135] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0136] 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.

[0137] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit it. Although the present invention has been described in detail with reference to the above embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the specific implementation of the present invention. Any modifications or equivalent substitutions that do not depart from the spirit and scope of the present invention should be covered within the scope of protection of the claims of the present invention.

Claims

1. A method for measuring the anti-turbulence heading stability of marine GNSS based on cross-epoch difference, characterized in that, Includes the following steps: S1. Real-time acquisition of raw observation data from shipborne GNSS receivers to construct a two-layer sliding buffer queue. The two-layer sliding buffer queue includes a long-interval buffer area for caching observation data of multiple historical epochs, and a short-interval buffer area for caching observation data of the current epoch and the previous adjacent epoch. S2, perform carrier phase difference between the current epoch and the previous adjacent epoch in the short interval buffer to calculate the short interval reference heading; S3. Select historical epochs that meet at least one of the preset time threshold and preset horizontal distance threshold from the long interval buffer as reference epochs. After performing ephemeris orbital error compensation and atmospheric delay compensation on the current epoch and the reference epoch, perform cross-epoch carrier phase differential to calculate the long interval coarse heading. S4, using the short-interval reference heading as the verification benchmark, perform heading autonomous integrity detection on the long-interval coarse heading to obtain the stable heading output for the current epoch; The autonomous heading integrity detection includes: calculating the heading difference between the short-interval reference heading and the long-interval coarse heading; if the heading difference exceeds a preset integrity threshold, the long-interval coarse heading is determined to be abnormal and the short-interval reference heading is output; otherwise, the long-interval coarse heading is output.

2. The method for measuring marine GNSS heading stability based on cross-epoch difference as described in claim 1, characterized in that, The calculation of the short-interval reference heading in step S2 includes: S21, Perform epoch-to-epoch difference on the carrier phase observations of the current epoch and the previous adjacent epoch to eliminate carrier integer ambiguity and solve for the relative position increment in the short-interval geocentric geofixed coordinate system; S22, the relative position increment in the short-interval geocentric-ground-fixed coordinate system is transformed to the ship's local northeast-sky coordinate system, the vertical undulation displacement is separated, the horizontal component is retained, and the short-interval eastward displacement increment and short-interval northward displacement increment are obtained. S23, perform heading calculation and quadrant correction on the short-interval eastward displacement increment and short-interval northward displacement increment to obtain the short-interval reference heading, and store it in the short-interval buffer as a verification benchmark for heading integrity detection.

3. The method for measuring marine GNSS heading stability based on cross-epoch difference as described in claim 2, characterized in that, Step S22 involves converting the relative position increment in the short-interval geocentric-fixed coordinate system to the ship's local northeast-sky coordinate system, separating the vertical undulation displacement, and retaining the horizontal component, including: S221, obtain the ship's latitude and longitude coordinates at the current epoch, and construct the rotation matrix from the geocentric earth-fixed coordinate system to the ship's local northeast-sky coordinate system; S222, using the rotation matrix, the relative position increment in the short-interval geocentric geofixed coordinate system is transformed to the ship's local northeast-sky coordinate system to obtain the eastward displacement component, the northward displacement component, and the skyward displacement component. S223, remove the celestial displacement component, and use the eastward displacement component and the northward displacement component as the short-interval eastward displacement increment and the short-interval northward displacement increment, respectively.

4. The method for measuring marine GNSS heading stability based on cross-epoch difference as described in claim 1, characterized in that, Step S3, after performing ephemeris orbital error compensation and atmospheric delay compensation on the current epoch and the reference epoch, performs cross-epoch carrier phase differential calculation to obtain the long-interval coarse heading, including: S321, Read the satellite broadcast ephemeris corresponding to the current epoch and the reference epoch respectively, calculate the satellite orbit coordinates and satellite clock difference at the two epochs, and calculate the line-of-sight orbit residual error; S322 uses the ionospheric model and the tropospheric model to calculate the atmospheric delay difference between two epochs; S323, using the line-of-sight orbital residual error and the atmospheric delay difference, pre-correct the carrier phase observation values ​​of the current epoch and the reference epoch; S324, perform cross-epoch difference on the pre-corrected carrier phase observations, solve for the long-interval relative position increment, and obtain the long-interval coarse heading through coordinate transformation and heading calculation.

5. The method for measuring marine GNSS heading stability based on cross-epoch difference as described in claim 4, characterized in that, Step S321, which involves calculating the satellite orbit coordinates and clock bias at two epochs and calculating the line-of-sight orbit residual error, includes: S3211, Read the satellite broadcast ephemeris corresponding to the current epoch and the reference epoch respectively, and calculate the satellite orbit coordinates and satellite clock difference at the two epochs; S3212, calculates the satellite orbital position deviation and satellite clock error difference between two epochs; S3213, Project the satellite orbital position deviation and the satellite clock error difference onto the satellite-receiver line-of-sight direction to obtain the line-of-sight orbital residual error.

6. The method for measuring marine GNSS heading stability based on cross-epoch difference as described in claim 4, characterized in that, Step S322, which uses an ionospheric model and a tropospheric model to calculate the atmospheric delay difference between two epochs, includes: S3221, Using the Klobuchar ionospheric model, the ionospheric delay between the current epoch and the reference epoch is calculated to obtain the ionospheric delay difference. S3222, using the Saastamoinen tropospheric model, calculate the tropospheric delay between the current epoch and the reference epoch to obtain the tropospheric delay difference; S3223, The ionospheric delay difference value and the tropospheric delay difference value are combined to obtain the atmospheric delay difference value.

7. The method for measuring marine GNSS heading stability based on cross-epoch difference as described in claim 4, characterized in that, Step S324, which involves performing cross-epoch differential analysis on the pre-corrected carrier phase observations to solve for the long-interval relative position increments and obtaining the long-interval coarse heading through coordinate transformation and heading calculation, includes: S3241 performs cross-epoch difference on the pre-corrected carrier phase observations to solve for the relative position increment in the long-interval geocentric-geofixed coordinate system. S3242, the relative position increment under the long-interval geocentric geofixed coordinate system is transformed to the ship's local northeast-sky coordinate system, and the vertical undulation displacement is separated to obtain the long-interval eastward displacement increment and the long-interval northward displacement increment. S3243, the long-interval coarse heading is obtained by calculating the heading of the long-interval eastward displacement increment and the long-interval northward displacement increment using the arctangent function.

8. The method for measuring marine GNSS heading stability based on cross-epoch difference as described in claim 1, characterized in that, Step S4, which uses the short-interval reference heading as a verification benchmark to perform heading autonomy integrity testing on the long-interval coarse heading to obtain a stable heading output for the current epoch, includes: S41, extract the short-interval reference heading from the short-interval buffer, obtain the long-interval coarse heading, and calculate the absolute value of the heading difference between the two; S42, compare the absolute value of the heading difference with the preset integrity threshold; if the absolute value of the heading difference is less than or equal to the preset integrity threshold, then take the long-interval coarse heading as the effective output heading; S43, if the absolute value of the heading difference is greater than the preset integrity threshold, it is determined that the long-interval coarse heading is abnormal, and the short-interval reference heading is output as the alternative heading for the current epoch. S44. Perform quadrant correction and normalization on the output heading determined in step S42 or step S43 to obtain stable heading data for the current epoch.

9. The method for measuring marine GNSS heading stability based on cross-epoch difference as described in any one of claims 1-8, characterized in that, After obtaining the stable heading output for the current epoch in step S4, the following steps are also included: S51, obtain the speed and heading values ​​for a consecutive preset number of observation epochs, and calculate the speed fluctuation value and heading change rate; S52, Based on the comparison results of the speed fluctuation value and the heading change rate with the preset speed fluctuation threshold and the preset heading change rate threshold, determine the current navigation condition of the ship. S53 adaptively adjusts the differential strategy and detection threshold based on the determined navigation conditions.

10. The method for measuring marine GNSS heading stability based on cross-epoch difference as described in claim 9, characterized in that, Step S53, which involves adaptively adjusting the differential strategy and detection threshold based on the determined navigation conditions, includes: S531, if it is determined that the ship is in a uniform speed straight navigation condition, continue to execute steps S2 to S4 in a loop, and use the long-interval coarse heading as the main output heading. S532, if it is determined that the ship has entered a turning or speed change maneuvering condition, the preset time threshold, the preset horizontal distance threshold and the preset integrity threshold are increased, and the maximum time difference of the reference epochs that can be screened in the long interval buffer is shortened; S533, under the steering or speed change maneuvering condition, reduce the output weight of the long-interval coarse heading and increase the output weight of the short-interval reference heading.

Images