Mid-orbit satellite orbit measurement data preprocessing method and device and computer equipment

Abstract

The application relates to a kind of middle orbit satellite orbit measurement data preprocessing method, device and computer equipment. By adopting layered progressive data cleaning mode, each stage can be independently executed test optimization, and through three-dimensional parallel deduplication and minimum uncertainty record reservation, the deduplication accuracy is significantly improved. In the aspect of outlier detection, three statistical methods are combined with physical constraint verification and time continuity check weighted voting, so that the false detection rate and the missed detection rate are significantly reduced. In the aspect of segmented discontinuity detection, the natural arc segment interval and the abnormal breakpoint are effectively distinguished, the false interpolation of the natural arc segment interval is avoided, and the physical meaning of the data is maintained. In the aspect of adaptive algorithm selection, the optimal filtering and resampling algorithm are automatically selected based on multidimensional feature analysis and the parameters are optimized layer by layer, without manual parameter adjustment, overcoming the limitations of fixed algorithm. Through four-dimensional comprehensive score, the data quality is comprehensively evaluated, and the data quality qualified rate is significantly improved.

Description

Technical Field

[0001] This application relates to the field of aerospace data processing technology, and in particular to a method, apparatus, computer equipment, computer-readable storage medium, and computer program product for preprocessing orbit measurement data of medium-Earth orbit satellites. Background Technology

[0002] In aerospace engineering, the preprocessing of orbit measurement data for Medium Earth Orbit (MEO) satellites is a crucial step in applications such as precise orbit determination, satellite constellation management, and collision warning. MEO relay satellites generally employ an X-band (standard frequency range of 8-12 GHz) ground station-transponder two-way pseudo-code ranging / velocity measurement system. This system calculates distance by measuring the round-trip time of the signal. This system can effectively eliminate ground station clock errors, but the measurement data is still affected by error sources such as tropospheric delay, transponder delay uncertainty, multipath effects, equipment noise, and ionospheric delay.

[0003] Raw measurement data often suffers from various quality issues during acquisition and transmission, including out-of-order timestamps, duplicate records, outliers, missing data, and segment discontinuities (including natural arc intervals caused by measurement geometry and abnormal breakpoints caused by non-geometric factors). If these issues are not addressed promptly, they will severely impact the accuracy and reliability of subsequent precise trajectory determination. Therefore, it is essential to systematically clean, filter, resample, and perform quality control on the raw measurement data.

[0004] In traditional technologies, the preprocessing methods widely used in aerospace engineering are mainly based on a fixed process: first, format verification and integrity checks are performed; then, outlier detection is performed using the 3σ criterion; data smoothing is performed using a simple moving average filter; finally, linear interpolation resampling is performed; and quality assessment is based on RMS (Root Mean Square) error. However, this type of method has significant shortcomings: First, it uses fixed algorithms and parameters, making it impossible to adaptively select the optimal algorithm based on data characteristics (such as signal-to-noise ratio, data distribution characteristics, sampling rate, etc.), resulting in processing effectiveness that relies on human experience and is inefficient. Second, the outlier detection method is singular; the 3σ criterion assumes that the data follows a normal distribution, which is not robust enough to non-normally distributed data and is prone to false positives or false negatives. Third, it cannot effectively identify and process the discontinuities between natural arc intervals and outlier segments in MEO satellite measurement data, which may lead to incorrect interpolation of natural intervals and damage to the physical meaning of the data. Fourth, the quality assessment system is incomplete, usually using only a single indicator such as RMS error, lacking a comprehensive assessment from multiple dimensions such as completeness, accuracy, consistency, and timeliness. Fifth, the system is mostly implemented using hard coding, lacking a plug-in architecture and configuration-driven mechanism, making it difficult to flexibly extend to new algorithms.

[0005] In addition, MEO satellite measurement data is characterized by large data volume, high real-time requirements, complex data characteristics, and strict quality requirements. Traditional tools and methods usually cannot meet the above engineering needs at the same time. Summary of the Invention

[0006] Therefore, it is necessary to provide a method, apparatus, computer equipment, computer-readable storage medium, and computer program product for preprocessing orbit measurement data of medium-Earth orbit satellites to address the aforementioned technical problems.

[0007] In a first aspect, this application provides a method for preprocessing orbit measurement data of a medium Earth orbit satellite, the method comprising: The acquired raw orbit measurement data is subjected to format recognition, integrity check and timestamp verification to obtain verified measurement data; the raw orbit measurement data includes at least one of X-band two-way ranging data, velocity data and GNSS auxiliary data from medium Earth orbit relay satellites; The validated measurement data is cleaned to obtain cleaned measurement data. The data cleaning process includes sequentially sorting the measurement data according to time series, detecting duplicate data, detecting outliers, and handling missing values. The boundary of the measurement arc segment is identified based on the time interval threshold. The cleaned measurement data is divided into at least one measurement arc segment. The theoretical visibility of the measurement arc segment is calculated in combination with the ground station location and orbit parameters to verify the division result of the measurement arc segment. This is to distinguish between the natural arc segment interval between adjacent measurement arc segments and the abnormal segment discontinuity within the measurement arc segment, and to take corresponding processing strategies. Multidimensional feature analysis in the time domain, frequency domain, and statistical domain is performed on the processed measurement data to obtain data features. Based on the data features, filtering algorithms and resampling algorithms are determined to filter, denoise, and resample the processed measurement data. The quality scores of the resampled measurement data are calculated from four dimensions: completeness, accuracy, consistency, and timeliness. The quality scores of each dimension are weighted and combined to obtain a comprehensive quality score. The data is then graded and processed according to the comprehensive quality score, and the preprocessed data is output.

[0008] In one embodiment, the step of performing format recognition, integrity check, and timestamp verification on the acquired raw orbit measurement data to obtain verified measurement data includes: The file header magic number of the original orbit measurement data is detected to identify the data format and version, and the parser corresponding to the data format and version is called to parse the original orbit measurement data; The parsed orbit measurement data is subjected to an integrity check, which includes three levels of checks: format level, semantic level, and reference level. The timestamps of the orbital measurement data that have passed the integrity check are unified to the unified time system of GPS or UTC. Leap second correction is performed based on the cumulative leap second difference between GPS and UTC, and the time interval between adjacent epochs is checked to ensure that it is consistent with the nominal sampling interval in order to obtain the verified measurement data.

[0009] In one embodiment, the duplicate data detection includes: The measurement data sorted according to time series is subjected to parallel duplicate detection from three dimensions: time, space, and attribute. The time dimension is used to detect whether the time difference between any two data records is less than a preset time threshold; the space dimension is used to detect whether the spatial distance between the two data records in a unified coordinate system is less than a preset distance threshold; and the attribute dimension is used to detect whether the difference between the pseudocode ranging observation values ​​of the two data records is less than a preset attribute threshold. The detection results from the three dimensions are aggregated through weighted voting to determine the duplicate data group, and the data record with the minimum uncertainty is retained within the duplicate data group.

[0010] In one embodiment, the outlier detection includes: Three statistical methods—Z-score, median absolute deviation, and interquartile range—were used to perform parallel anomaly detection on the observation sequence of the measurement data after duplicate data detection, in order to determine the anomaly detection result for each sequence point. Physical constraints are verified for each sequence point based on the physical constraints of medium Earth orbit satellite orbit dynamics to determine the anomaly verification result for each sequence point; the physical constraints include position constraints, velocity constraints, and acceleration constraints. The observation sequence is checked for temporal continuity to determine the temporal continuity result for each sequence point. The anomaly detection result, anomaly verification result, and temporal continuity result of each sequence point are weighted and voted on according to preset weights. When the weighted score exceeds the preset score threshold, the corresponding sequence point is determined to be an anomaly.

[0011] In one embodiment, the missing value processing includes: A theoretical epoch sequence is generated based on the start and end times of the measurement data and the nominal sampling interval. The theoretical epoch sequence is matched with the measured epochs one by one, and the unmatched theoretical epochs are identified as missing data. When the length of the missing interval is less than the first preset number of epochs, cubic spline interpolation is used to interpolate the missing interval to obtain the interpolated epoch sequence. When the length of the missing interval is not less than the first preset epoch number and not greater than the second preset epoch number, the missing interval is interpolated by linear interpolation or orbit extrapolation to obtain the interpolated epoch sequence. When the length of the missing interval is greater than the second preset epoch number, the missing interval is marked as an unusable segment and no interpolation is performed.

[0012] In one embodiment, the detection of abnormal segment discontinuities employs a multi-method fusion strategy combining gradient detection, Pettitt test, and Mann-Kendall test; the identified abnormal segment discontinuities are processed independently, smoothly transitioned, or marked and preserved; and the identified natural arc intervals are left uninterpolated.

[0013] In one embodiment, the step of performing multidimensional feature analysis on the processed measurement data in the time domain, frequency domain, and statistical domain to obtain data features, and determining a filtering algorithm based on the data features, includes: The noise level is estimated based on the energy of the high-frequency components of the measurement data, and the signal-to-noise ratio is calculated; a normalized smoothness metric is calculated based on the second-order difference of the observation sequence; periodicity detection and non-stationarity detection are performed, and the non-stationarity detection determines whether the measurement data is stationary by calculating the rate of change of the variance of the sliding window. When the signal-to-noise ratio is greater than the preset signal-to-noise ratio threshold and the smoothness metric is less than the preset smoothness threshold, the filtering algorithm is determined to be Savitzky-Golay filtering. When the signal-to-noise ratio is not greater than the preset signal-to-noise ratio threshold, the filtering algorithm is determined to be moving average filtering; when the data has outlier characteristics, the filtering algorithm is determined to be median filtering. When the proportion of high-frequency components is greater than the preset proportion threshold, the filtering algorithm is determined to be low-pass filtering.

[0014] In one embodiment, the weighted combination of the quality scores of each dimension to obtain the comprehensive quality score includes: Where C represents the completeness score and A represents the accuracy score. For consistency scoring, T represents timeliness scoring. , , and For the corresponding weights, This is the overall quality score.

[0015] Secondly, this application also provides a medium Earth orbit satellite orbit measurement data preprocessing device, the device comprising: The data access module is used to perform format recognition, integrity check, and timestamp verification on the acquired raw orbit measurement data to obtain verified measurement data; the raw orbit measurement data includes at least one of X-band two-way ranging data, velocity data, and GNSS auxiliary data from medium Earth orbit relay satellites; The data cleaning module is used to clean the validated measurement data to obtain cleaned measurement data. The data cleaning process includes sequentially performing the measurement data sorting according to time series, duplicate data detection, outlier detection, and missing value handling. The arc segment identification module is used to identify the boundary of the measurement arc segment based on the time interval threshold, divide the cleaned measurement data into at least one measurement arc segment, and verify the division result of the measurement arc segment by combining the ground station position and orbit parameters to calculate the theoretical visibility, so as to distinguish the natural arc segment interval between adjacent measurement arc segments and the abnormal segment discontinuity within the measurement arc segment, and take corresponding processing strategies. The data filtering module is used to perform multi-dimensional feature analysis of the processed measurement data in the time domain, frequency domain, and statistical domain to obtain data features. Based on the data features, the filtering algorithm and resampling algorithm are determined to filter, denoise, and resample the processed measurement data. The data output module is used to calculate the quality score of the resampled measurement data from four dimensions: completeness, accuracy, consistency and timeliness. The quality scores of each dimension are weighted and combined to obtain the comprehensive quality score. The data is then graded and processed according to the comprehensive quality score, and the preprocessed data is output.

[0016] Thirdly, this application also provides a computer device, including a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program to implement the steps of the method described in the first aspect.

[0017] Fourthly, this application also provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the steps of the method described in the first aspect.

[0018] Fifthly, this application also provides a computer program product, including a computer program that, when executed by a processor, implements the steps of the method described in the first aspect above.

[0019] The aforementioned preprocessing methods, devices, computer equipment, computer-readable storage media, and computer program products for medium-Earth orbit satellite orbit measurement data employ a hierarchical, progressive data cleaning pipeline model. Each stage can be independently tested and optimized, and the accuracy of deduplication is significantly improved through three-dimensional parallel deduplication and retention of records with minimum uncertainty. In outlier detection, three statistical methods are used in parallel, combined with physical constraint verification and time continuity checks using weighted voting, significantly reducing both false positive and false negative rates, effectively overcoming the shortcomings of the traditional single 3σ criterion. In segmented discontinuity detection, natural arc intervals and abnormal breakpoints are effectively distinguished, avoiding erroneous interpolation of natural arc intervals and preserving the physical meaning of the data. Regarding adaptive algorithm selection, the optimal filtering and resampling algorithms are automatically selected based on multi-dimensional feature analysis, and parameters are optimized hierarchically, eliminating the need for manual parameter tuning and overcoming the limitations of fixed algorithms. A comprehensive four-dimensional scoring system is used to evaluate data quality, increasing the data quality pass rate to over 95%. Furthermore, thanks to the layered and progressive processing flow and parallel verification mechanism, the processing time is significantly shortened. Under a 32GB RAM and Intel Xeon (8 cores) environment, the processing latency for 10k points is <1s (recommended <100ms), meeting the requirements for near real-time processing and providing a high-quality data foundation for subsequent applications such as precise orbit determination, satellite constellation management, and collision warning. Attached Figure Description

[0020] To more clearly illustrate the technical solutions in the embodiments of this application or related technologies, the drawings used in the description of the embodiments of this application or related technologies will be briefly introduced below. Obviously, the drawings described below are only some embodiments of this application. For those skilled in the art, other related drawings can be obtained based on these drawings without creative effort.

[0021] Figure 1 This is a flowchart illustrating a method for preprocessing orbit measurement data of a medium Earth orbit satellite in one embodiment; Figure 2 This is a flowchart illustrating the quality score weighting process in one embodiment; Figure 3 This is a structural block diagram of a medium-Earth orbit satellite orbit measurement data preprocessing device in one embodiment; Figure 4 This is an internal structural diagram of a computer device in one embodiment. Detailed Implementation

[0022] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description is provided in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the scope of this application.

[0023] It should be noted that the terms "first," "second," etc., used in this application can be used to describe various elements, but these elements are not limited by these terms. These terms are only used to distinguish the first element from the second element. The terms "comprising" and "having," and any variations thereof, used in this application, are intended to cover non-exclusive inclusion. The term "multiple" used in this application refers to two or more. The term "and / or" used in this application refers to one of the embodiments, or any combination of multiple embodiments.

[0024] In one exemplary embodiment, this application provides a preprocessing method for medium Earth orbit satellite orbit measurement data, applicable to near real-time preprocessing of X-band two-way ranging data (ground station transmits signals, which are relayed by the satellite transponder and returned; the round-trip time is measured to calculate the distance) / velocity data of medium Earth orbit relay satellites, as well as GNSS (Global Navigation Satellite System) auxiliary data. This method employs hierarchical progressive data cleaning, segmented discontinuity detection, adaptive algorithm selection, and multi-dimensional quality assessment, and achieves flexible, configuration-driven processing through a plug-in architecture.

[0025] like Figure 1 As shown, the method may specifically include the following steps: Step 102: Perform format recognition, integrity check and timestamp verification on the acquired raw orbit measurement data to obtain verified measurement data.

[0026] The raw orbit measurement data includes at least one of the following: X-band two-way ranging data, velocity data, and GNSS auxiliary data from a medium Earth orbit relay satellite. This embodiment uses X-band two-way pseudo-code ranging data from a MEO relay satellite (orbital altitude h = 11500 km) as an example to illustrate its data preprocessing flow. This satellite employs an X-band ground station-transponder two-way pseudo-code ranging / velocity measurement system, with an uplink frequency of 7.2 GHz and a downlink frequency of 8.4 GHz, a frequency ratio of 8.4 / 7.2 = 7 / 6, which meets the standard transponder ratio. Its two-way system effectively eliminates ground station clock errors and improves ranging accuracy.

[0027] For example, the raw orbit measurement data in this embodiment is X-band two-way ranging data, using RINEX 3.05 format. The file consists of a header and epoch observation records. The header includes fields such as version and type, station / antenna / receiver information, approximate geocentric coordinates, system and observation type list, first and last observation times, and time system. The epoch record provides the time, epoch marker, number of satellites in this epoch, and lists the two-way pseudocode distance ρ (m) and Doppler / distance change rate observations by satellite. The data time range is from 2025-12-15T00:00:00 to 2025-12-15T02:46:40 (UTC), the nominal sampling rate is 1Hz, and the theoretical number of epochs N_theory=T / Δt_s=10000, where T represents the duration of the time period, which is 10000s in this embodiment; Δt_s represents the nominal sampling interval, which is 1s in this embodiment.

[0028] In this embodiment, after acquiring the raw orbit measurement data, format identification is performed first, for example, by detecting the magic number in the file header to identify the data format and version. If this embodiment identifies it as RINEX 3.05 format and version, then the parser corresponding to the RINEX 3.05 format and version is invoked to parse the data. Here, RINEX stands for Receiver Independent Exchange Format.

[0029] Following this, data integrity checks and timestamp verification are performed. The integrity check employs a three-level mechanism: format-level, semantic-level, and reference-level. The format-level check examines required header fields, including version, type, observation type table, and start and end times. The semantic-level check verifies epoch-by-epoch whether the number of fields matches the header definition, whether the satellite system / number is valid, whether the timestamp is resolvable and falls within a reasonable range, and whether the pseudorange and velocity values ​​are not NaN and whether their magnitudes conform to MEO X-band ranging characteristics (e.g., round-trip distance approximately 2 × 10⁻⁶). 4 (On the order of km); reference-level checks further validate the data using reference information. During integrity checks, processing terminates and diagnostic information is output when a fatal error exists in the header segment; single-epoch misalignments are skipped and logged. The data missing rate is defined as R_miss = (N_theory) N_obs) / N_theory×100%, in this example, the time period is T=10000 s, the nominal sampling interval is Δt_s=1 s, and the theoretical number of epochs is N_theory=T / Δt_s=10000. After parsing the original record N_raw=9758, and removing 8 duplicates, the result is still 9750, therefore N_obs=9750. After aligning the grid, the missing data N_miss=250, so the missing rate R_miss=2.5%. If the missing rate threshold is 5%, since the missing rate R_miss=2.5% is less than the threshold of 5%, it means the integrity check has passed.

[0030] For timestamp verification, the timestamps of data that pass the above integrity checks are unified to GPS time or UTC unified time system; leap second correction is performed according to t_UTC=t_GPS. Δt_LS is executed, where Δt_LS is the accumulated GPS-UTC leap second difference since January 6, 1980. Since January 1, 2017, Δt_LS = 18 s. In this embodiment, the data date is December 15, 2025, so Δt_LS = 18 s is taken. The GPS timestamp is converted to UTC and compared with the nominal timeline, with a tolerance based on the receiver resolution (milliseconds). Simultaneously, the time interval t_i+1 between adjacent epochs after sorting is verified. The consistency between t_i and the nominal sampling interval allows for integer multiple intervals due to missing data (in this embodiment, the timestamp is strictly monotonically increasing), ultimately resulting in verified measurement data.

[0031] Step 104: Perform data cleaning on the verified measurement data to obtain cleaned measurement data.

[0032] The data cleaning process includes sequentially sorting the measurement data according to time series, detecting duplicate data, detecting outliers, and handling missing values.

[0033] For example, time series sorting is performed first, which can be done using the Timsort algorithm. For instance, sorting keys can include primary key (timestamp), secondary key (satellite ID), and tertiary key (observation station ID). After sorting, the data is strictly arranged in ascending order of timestamps, with a time complexity of O(n log n), where n = 10000 is the number of data records. This sorting lays the foundation for subsequent time series-based cleaning operations.

[0034] Then, duplicate data detection is performed, checking any two data records (not necessarily adjacent) in parallel across three dimensions: time, space, and attribute. Specifically, the time dimension detection checks whether the time difference between any two data records in the dataset is less than a preset time threshold, criterion being |t_i|. t_j|<10 -6s, where t_i and t_j are the timestamps of data i and j in any two data sets. For example, in this embodiment, a hash table is built based on the quantized time key, and a full table query is performed on the above data, detecting 5 duplicate data points.

[0035] Spatial dimension detection is performed by calculating whether the spatial distance between two data records in a unified coordinate system (ECEF or ECI) is less than a preset distance threshold, with the criterion being d = ||r_i||. The threshold r_j‖2 < 0.1 m is used to capture repeated sampling of nearly identical geometric points, where r_i and r_j are the three-dimensional positions of data i and data j in any two data sets. Pure two-way ranging files typically do not contain instantaneous geocentric positions; r is obtained from coarse ephemeris, geometric calculations, or GNSS-assisted states. If r has not yet been generated, this dimension detection can be skipped or downgraded and noted in the metadata. For example, this embodiment detects two duplicate data points based on this.

[0036] Attribute dimension detection checks whether the difference between the pseudocode ranging observations of two data records is less than a preset attribute threshold, with the criterion being |ρ_i ρ_j|<ε_ρ (precise repeatability up to 10) -6 m (approximate repetition can be relaxed according to ranging noise, such as 0.5m, subject to configuration). Here, ρ_i and ρ_j are the pseudo-code ranging values ​​of data i and data j from any two data sets. It should be noted that the pseudo-code ranging value ρ is a distance-type observation obtained by multiplying the pseudo-random code phase delay by the speed of light. For example, this embodiment detects 3 duplicate data points based on this.

[0037] After parallel detection of duplicate data across three dimensions—time, space, and attributes—the results of the duplicate detection across these three dimensions are aggregated using a weighted voting strategy to determine the duplicate data group. Within the duplicate data group, the data record with the lowest uncertainty is retained, i.e., i*=arg min σ_i. The uncertainty is preferentially determined using the original built-in precision information (RINEX sidebar SSI / LLI, or σ_ρ, quality term provided by the device), and only when the original information is missing is it calculated using the SNR-σ empirical model or error propagation. For example, in this embodiment, 8 duplicate data points are ultimately removed, leaving 9750 valid records.

[0038] Outlier detection is then performed, specifically using three statistical methods: Z-score, Median Absolute Deviation (MAD), and Interquartile Range (IQR). These methods are used to perform parallel anomaly detection on the observation sequence of the measurement data after duplicate data detection. The observation sequence x can be a pseudorange, velocity measurement, or derived kinematic components. For example, the Z-score (standard score, the deviation of a data point from the mean normalized to the standard deviation) method calculates Z = (x...)... μ) / σ, where μ and σ are the sample mean and standard deviation, and |Z|>3 is a candidate anomaly (normal distribution covers approximately 99.7%). This embodiment detected 12 candidate anomaly sequence points.

[0039] The MAD method calculates MAD = median(|x_i) by... median(x)|), correct Z_MAD=0.6745·(x) The median(x) / MAD criterion |Z_MAD|>3 is used to identify candidate anomalies, and the coefficient 0.6745 makes Z_MAD under normal distribution comparable to the ordinary Z-score scale. This embodiment detected 15 candidate anomaly sequence points.

[0040] The IQR method involves taking the first quartile (Q1) and the third quartile (Q3) and calculating IQR = Q3. Q1, falling on [Q1] Any value outside of [1.5·IQR, Q3+1.5·IQR] is considered an anomaly (Tukey box rule). This embodiment detected 10 candidate anomaly sequence points.

[0041] After outlier detection using three statistical methods—Z-score, median absolute deviation (MAD), and interquartile range (IQR)—anomaly sequence points are determined through a voting mechanism. For example, sequence points that are detected as candidate anomalies by at least two of the three methods in parallel detection are identified as anomaly sequence points. In this embodiment, 18 anomaly sequence points were obtained.

[0042] It is understandable that, in addition to the three statistical methods listed above—Z-score, MAD, and IQR—sliding window detection, machine learning-based anomaly detection, and wavelet transform-based anomaly detection methods can also be used. For specific application scenarios, appropriate anomaly detection methods or combinations thereof can be selected based on data characteristics and accuracy requirements; this embodiment does not impose any limitations on this.

[0043] Further physical constraint verification was performed on all sequence points. Specifically, physical constraint verification for each sequence point was performed based on the physical constraints of the orbital dynamics of medium Earth orbit satellites (scanning all data points, not just checking marked outliers). These physical constraints include: position constraints, geocentric distance r ∈ [11000, 14500] km; velocity constraints, velocity modulus ‖v‖ ∈ [4.5, 5.5] km / s (typically the theoretical velocity is approximately 4.72 km / s in a circular orbit with h = 11500 km); and acceleration constraints, non-conservative force acceleration < 10. -6 km / s². It should be noted that the velocity is preferably obtained from the original Doppler; otherwise, it is obtained from the position difference.

[0044] Simultaneously, a time continuity check is performed on the observed sequence based on temporal continuity, such as checking whether the rate of change between adjacent epochs is within a reasonable range. Finally, the statistical anomaly detection results, physical constraint verification results, and time continuity check results for each sequence point are weighted and voted on according to preset weights. The weights for statistical anomaly detection are 0.4, physical constraint check are 0.4, and time continuity check are 0.2. When the weighted score exceeds a preset score threshold (e.g., 0.5), the corresponding sequence point is determined as the final anomaly. In this embodiment, 18 anomaly sequence points (i.e., statistical anomaly detection results) are identified through parallel detection using the three statistical methods. Based on the physical constraint verification results and time continuity check results, a weighted vote is performed according to preset weights, and sequence points with a weighted score greater than 0.5 are determined as the final anomaly sequence points. In this embodiment, 15 final anomaly sequence points are obtained.

[0045] Finally, missing values ​​are handled. Specifically, a theoretical epoch sequence is generated based on the start and end times of the measured data and the nominal sampling interval. An expected epoch sequence is generated from t_0, t_1, and Δt_s. The theoretical epoch sequence is then matched one by one with the measured epochs. If |t_obs If t_exp|<τ_tol (usually taken as 0.5·Δt_s, 0.5 s in this embodiment), then a match is found; otherwise, a mismatch is recorded and a missing data point is identified, thus determining the unmatched theoretical epoch as missing data. Based on this, this embodiment detected 250 missing data points, with a missing rate of 2.5%. When the length of the missing interval is less than the first preset epoch number (e.g., 3 epochs), cubic spline interpolation is used for the missing interval; in this embodiment, this results in 180 points. When the length of the missing interval is not less than (i.e., greater than or equal to) the first preset epoch number and not greater than (i.e., less than or equal to) the second preset epoch number (e.g., 10 epochs), linear interpolation or orbital extrapolation is used for the missing interval; in this embodiment, this results in 60 points. When the length of the missing interval is greater than the second preset epoch number (e.g., 10 epochs), the missing interval is marked as an unusable segment and no interpolation is performed; in this embodiment, this results in 10 points.

[0046] Calculate the uncertainty after interpolation Where σ_raw is the nearest measured uncertainty, Δt is the interval (s) between the interpolation time and the most recent measured time, and k is the model error growth coefficient (s). - ¹ (calibrated by cross-validation). This formula is used to amplify the uncertainty of the interpolation point with Δt; a quality label "interpolated" is added to the interpolated data to distinguish it from the measured data, ensuring accurate perception of data quality in subsequent processing. This completes the data cleaning process.

[0047] This embodiment employs three statistical methods—Z-score, MAD, and IQR—in parallel detection, combined with weighted voting confirmation based on physical constraints of MEO satellite orbital dynamics (geocenter distance, velocity modulus, and acceleration) and time continuity checks. Compared to traditional single-method approaches that rely solely on the 3σ criterion, this embodiment demonstrates stronger robustness to non-normally distributed data. Its multi-strategy fusion mechanism significantly reduces both the false positive and false negative rates of outlier detection.

[0048] Step 106: Identify the boundary of the measurement arc segment based on the time interval threshold, divide the cleaned measurement data into at least one measurement arc segment, and perform discontinuity detection and processing.

[0049] Here, an arc segment represents the measurement period from satellite rise to fall (Loss of Signal), with the visible arc segment length for a single MEO satellite station being approximately 30-90 minutes. Segment discontinuity refers to significant segment discontinuities that occur during data acquisition and transmission, including natural arc segment intervals (caused by measurement geometry) and anomalous segment discontinuities (caused by non-geometric factors).

[0050] In this embodiment, the boundary of the measurement arc segment is identified based on the time interval threshold. The cleaned measurement data is divided into at least one measurement arc segment. The theoretical visibility of the measurement arc segment division is verified by combining the ground station location and orbit parameters to distinguish between the natural arc segment interval between adjacent measurement arc segments and the abnormal segment discontinuity within the measurement arc segment, and corresponding processing strategies are adopted.

[0051] For example, the boundaries of measurement arcs are first automatically identified based on a time interval threshold. In this embodiment, a time interval threshold θ = 600 s (10 min) is used as an example. When the time difference between adjacent existing epochs exceeds this threshold, a new measurement arc is cut off. The effective arc length is usually determined by the actual duration of single-station visible tracking (such as about 30-90 min for MEO satellites), and is constrained by measurement geometry, so it does not need to be close to the threshold θ. Based on this, the cleaned measurement data is divided into at least one measurement arc, and the theoretical visibility is calculated in combination with the ground station position and orbital parameters. The division results of the measurement arc are verified to accurately distinguish between the natural arc intervals between adjacent measurement arcs (caused by measurement geometry) and the abnormal segmentation discontinuities within the measurement arc (caused by non-geometric factors), and corresponding processing strategies are adopted.

[0052] In this embodiment, three measurement arcs were detected: arc 1 (00:00:00-00:45:00, approximately 2700 data points), arc 2 (01:30:00-02:15:00, approximately 2700 data points), and arc 3 (02:30:00-02:46:40, approximately 1000 data points). Adjacent data gaps were 45 min (2700 s) and 15 min (approximately 900 s), both greater than 600 s, thus allowing for stable segmentation into three segments. The arc segment division results were verified by combining the ground station location and orbital parameter calculations to ensure theoretical visibility. This confirmed that the aforementioned gaps were natural arc intervals (caused by measurement geometry). Therefore, during processing, the intervals were not interpolated; only the arc boundaries were marked and passed to the subsequent orbital determination algorithm, thereby avoiding erroneous interpolation that could corrupt the physical meaning of the data.

[0053] Then, anomaly segmentation detection is performed. This embodiment employs a multi-method fusion strategy combining gradient detection, Pettitt's test, and Mann-Kendall's test. An anomaly segment is only identified if at least two of these methods detect anomalies. Specifically, gradient detection is used for samples observed both inside and on both sides of the arc segment, calculating the first-order difference Δxi = x(i+1). The difference between xi and the second-order difference Δ²xi is used to obtain candidate points with |Δ²xi|>k·σ_Δ² (default k=3), and adjacent candidates that are too close are removed. The difference is not calculated for time windows with no observations. For example, an anomalous breakpoint 01:48:30 is detected within arc 2 (01:30:00-02:15:00), with an abnormally large gradient.

[0054] A physical rationality check is then performed to verify whether the jump variable level exceeds the task threshold, whether it is adjacent to track control / station switching / obstruction, and whether it is consistent with visibility abrupt changes. If it matches a known maneuver, it is marked as an event boundary; if it is statistically significant and has no physical correspondence, it is confirmed as an abnormal breakpoint. In this embodiment, the breakpoint at 01:48:30 was confirmed as an abnormal breakpoint after the physical rationality check. Furthermore, the identified discontinuous abnormal segments are handled using a segmented independent processing, smooth transition, or mark retention strategy. This embodiment employs a segmented independent processing strategy for abnormal breakpoints, cleaning and filtering the abnormal segments separately, while keeping the identified natural arc intervals uninterpolated, thus ensuring the rigor of data processing.

[0055] Step 108: Perform multidimensional feature analysis of the processed measurement data in the time domain, frequency domain, and statistical domain to obtain data features, and determine the filtering algorithm and resampling algorithm based on the data features.

[0056] In this embodiment, multidimensional feature analysis of the processed measurement data is performed in the time domain, frequency domain, and statistical domain to obtain data features. Based on the data features, filtering algorithms and resampling algorithms are determined to filter, denoise, and resample the processed measurement data.

[0057] For example, the processed measurement data is first subjected to multidimensional feature analysis in the time domain, frequency domain, and statistical domain to obtain data features. For instance, in terms of noise level assessment, the noise level is estimated based on the energy of the high-frequency components of the measurement data and the signal-to-noise ratio (SNR) is calculated. In this embodiment, the SNR is found to be 35 dB.

[0058] Regarding the smoothness metric for signal feature recognition, the normalized smoothness metric is calculated based on the second-order difference of the observed sequence, and can be specifically calculated using the following formula (1): (1) in =std(x), which is the standard deviation of the observed sequence; For a dimensionless parameter, the criterion is: <ε1, in this embodiment ε1=0.05, actual measurement =0.012<ε1, therefore it is judged as relatively smooth.

[0059] In terms of periodic detection, the orbital period was identified as approximately 6.6 hours.

[0060] In nonstationarity detection, the rate of change of variance of the sliding window, Δσ² / σ² = |Var(Wk+1), is calculated. The expression Var(Wk)| / Var(Wk) determines whether the measurement data is stationary, that is, whether the statistical characteristics of the measurement data remain stable over time. Generally, a rate of change less than 0.5 is considered stationary. Here, Var(Wk) and Var(Wk+1) represent the sample variances of the observations within the k-th and (k+1)-th sliding windows, respectively, reflecting the degree of dispersion of the data within that time window. In this example, Δσ² / σ² = 0.15 < 0.5, therefore it is judged to be stationary.

[0061] Then, the filtering algorithm is determined based on the above data characteristics. This embodiment uses a rule-based expert system / decision tree to select the filtering algorithm. For example, when the signal-to-noise ratio (SNR) is greater than a preset SNR threshold and the smoothing metric is less than a preset smoothing threshold, the filtering algorithm is determined to be Savitzky-Golay filtering. In this embodiment, SNR = 35 dB > 30 dB and... Since ε < 1, the Savitzky-Golay filter is chosen, with a recommended window size w = 11 and polynomial order n = 3.

[0062] In addition, if the signal-to-noise ratio is not greater than (i.e. less than or equal to) a preset threshold, a moving average filter is selected; if outlier features exist, a median filter is selected; if the proportion of high-frequency components is greater than a preset proportion threshold, a low-pass filter (Butterworth) is selected; and under strong real-time constraints, methods with smaller computational load are preferred.

[0063] Further parameter optimization can be performed, specifically using a hierarchical parameter optimization method. For example, a coarse optimization (grid search) is first performed, searching for the optimal combination in the parameter space {window=[9,11,13], poly=[2,3,4]}. In this embodiment, the optimal parameters obtained after grid search are window=11 and poly=3. Subsequently, fine optimization (such as Bayesian optimization) and validation optimization (such as K-fold cross-validation) are performed, without requiring manual parameter tuning throughout the process. Finally, a resampling algorithm (such as cubic spline interpolation, linear interpolation, nearest neighbor interpolation, piecewise cubic Hermite interpolation (PCHIP), B-spline interpolation, or interpolation based on the orbital dynamics model) is determined based on the data characteristics. The filtered and denoised measurement data is reconstructed onto a uniform, equally spaced standard epoch grid (such as 1 Hz sampling), thereby outputting a time-continuous and uniformly sampled observation sequence to meet the time regularity requirements of the subsequent orbit determination algorithm.

[0064] This embodiment determines the filtering algorithm based on data characteristics, rather than using a fixed filtering algorithm as in traditional techniques. This not only avoids the problem of over-smoothing or under-smoothing of data with moderate noise levels using fixed algorithms, but also effectively preserves the first and second derivative characteristics of the signal while eliminating high-frequency noise. Coarse parameter optimization through grid search ensures that the selected parameters are optimal in the local parameter space, improving processing efficiency by more than 30% compared to manual parameter tuning.

[0065] Step 110: Calculate the quality score of the resampled measurement data from four dimensions: data completeness, accuracy, consistency, and timeliness, and then classify the data according to the quality score.

[0066] In this embodiment, as Figure 2 As shown, the quality scores of the resampled measurement data are calculated from four dimensions: data integrity, accuracy, consistency, and timeliness. The quality scores of each dimension are weighted and combined to obtain the comprehensive quality score. The data is then graded and processed according to the comprehensive quality score, and the preprocessed data is output.

[0067] The integrity score is calculated as the data integrity rate C = 100 × N_valid / N_theory, where N_valid represents the number of valid data points and N_theory represents the theoretical epochs. In this embodiment, N_valid = 9990 and N_total = 10000, so substituting these values ​​into the calculation yields C = 100 × 9990 / 10000 = 99.90 points.

[0068] Accuracy scoring is determined by calculating the root mean square error relative to a reference value. A = 100 × (1 RMSE / Range), where the root mean square error is... In the formula, N represents the total number of epochs of measurement data involved in the evaluation. This represents the actual observed value at each K epoch, such as a pseudocode ranging value or a velocity measurement value. This represents the reference value corresponding to the k-th epoch, and is also used to compare it with the actual observed value. The benchmark value for comparison. Range is the range of the evaluated scalar (such as pseudocode ranging, ranging residual), Range = max(yk). min(yk) represents the number of participants. This is a distance-type or residual-type scalar, such as pseudocode ranging or ranging residual. In this embodiment, RMSE=0.15m and Range=500m, so substituting the parameters, we get A=99.97 points.

[0069] Consistency scoring is achieved by calculating the outlier ratio Co = 100 × (1 N_outlier / N_valid). Where N_outlier represents the number of outliers. In this embodiment, the number of outliers N_outlier=15 and N_valid=9990, so substituting the parameters, Co=100×(1-15 / 9990)=99.85 points.

[0070] Timeliness scoring is achieved by calculating the time-delay decay of data, T = 100 × exp( Δt / τ), where Δt=2 s is the delay from the measurement epoch to the preprocessing output, and τ=3600 s is the system design time scale, so T≈99.94 points.

[0071] The overall quality score is obtained by weighting and combining the quality scores of each dimension. For example, the following formula (2) can be used for weighting and combining: (2) in, The overall quality score is calculated using weighted averages, where C represents the completeness score and A represents the accuracy score. For consistency scoring, T represents timeliness scoring. , , and These are corresponding weights, specifically configurable values. Each weight can be dynamically adjusted based on the application scenario; for example, the timeliness weight can be increased for real-time processing scenarios. For historical analysis scenarios, timeliness can be ignored, such as configuration. =0.00, which can increase the accuracy weight for high-reliability applications. Consistency weight .

[0072] For example, in this embodiment, It is 0.25. It is 0.35. It is 0.25. Taking 0.15 as an example, the overall quality score is... .

[0073] Then, the data is graded based on the overall quality score. For example, a score of 90-100 is defined as excellent and can be directly used as the preprocessed data output; a score of 80-89 is defined as good, which is only applicable when at most one single item score is below 80. Local parameter adjustments are made to the selected algorithm and the cleaned results (such as fine-tuning the filter window length by ±2 to maintain an odd number, adjusting the Butterworth cutoff frequency by ±10%, replacing labeled isolated outliers with the neighborhood median, and fine-tuning interpolation parameters for short missing segments). The dimension with the lowest score is identified first, and then the relevant parameters for that dimension are adjusted. This process is repeated a maximum of two rounds. If the re-evaluation reaches a certain level... Output the results after the score is ≥90 and each individual score is not lower than 85; when the overall quality score is (70-79 points), it is defined as the general level. If it is determined to be at this level, feature analysis needs to be performed again and global parameter tuning is required for parameter optimization; when the overall quality score is (<70 points), it is defined as the poor level. If it is determined to be at this level, return to the preprocessing starting point and clean again.

[0074] In this embodiment =99.92, which is considered excellent, so the preprocessed data is directly output. This means that after the complete data preprocessing process, the original 10,000 theoretical epochs output 9,990 high-quality data points (including interpolation backfilling), with a comprehensive data quality score of 99.92, meeting the requirements for high-precision orbit determination. In a 32GB RAM, Intel Xeon (8 cores) environment, the total processing time is approximately 850ms, meeting near real-time processing requirements (<1s).

[0075] In the aforementioned preprocessing method for medium Earth orbit satellite orbit measurement data, a hierarchical, progressive data cleaning pipeline is adopted, allowing each stage to be independently tested and optimized. Three-dimensional parallel deduplication and retention of records with minimum uncertainty significantly improve deduplication accuracy. For outlier detection, three statistical methods are used in parallel, combined with physical constraint verification and time continuity checks using weighted voting, significantly reducing both false positive and false negative rates, effectively overcoming the shortcomings of the traditional single 3σ criterion. For segmented discontinuity detection, natural arc intervals and abnormal breakpoints are effectively distinguished, avoiding erroneous interpolation of natural arc intervals and preserving the physical meaning of the data. Regarding adaptive algorithm selection, the optimal filtering and resampling algorithms are automatically selected based on multi-dimensional feature analysis, and parameters are optimized hierarchically, eliminating the need for manual parameter tuning and overcoming the limitations of fixed algorithms. A comprehensive four-dimensional scoring system is used to fully evaluate data quality, increasing the data quality pass rate to over 95%. Furthermore, thanks to the layered and progressive processing flow and parallel verification mechanism, the processing time is significantly shortened. Under a 32GB RAM and Intel Xeon (8 cores) environment, the processing latency for 10k points is <1s (recommended <100ms), meeting the requirements for near real-time processing and providing a high-quality data foundation for subsequent applications such as precise orbit determination, satellite constellation management, and collision warning.

[0076] In one exemplary embodiment, the preprocessing of onboard GNSS-assisted data from a MEO relay satellite is used as an example. The complete preprocessing flow of this embodiment includes data input verification, data cleaning, segment discontinuity detection, adaptive algorithm selection, and multidimensional quality assessment.

[0077] In the data input verification phase, the input data consisted of time-series position and velocity data processed by the spaceborne GNSS receiver and transmitted via telemetry. The data was in RINEX 4.02 format, containing 5000 observation epochs, a sampling rate of 0.1 Hz, and a time range from 2025-12-15T00:00:00 to 2025-12-15T13:53:20 (UTC). The system identified the file as RINEX 4.02 format using the magic number in the file header and called the corresponding parser for parsing. After format verification, the data integrity check showed a missing rate of 1.2% (less than the 5% threshold), thus passing the check. Timestamp verification also passed.

[0078] During the data cleaning phase, the Timsort algorithm was used for time series sorting, and the sorted data was strictly arranged in ascending order of timestamps. Duplicate data detection was performed in parallel across three dimensions: time, space, and attributes. Two duplicate data points were found, and the strategy of retaining the record with the minimum uncertainty was used to remove duplicates. Outlier detection employed three statistical methods in parallel: Z-score, MAD, and IQR, combined with physical constraint checks (location constraint: geocentric distance r ∈ [11000, 14500] km, velocity constraint: velocity modulus ||v|| ∈ [4.5, 5.5] km / s), detecting a total of 8 outliers. An outlier correction strategy was adopted (e.g., replacing isolated outliers with the median of adjacent values). Missing value processing detected 60 missing data points. Short-term missing data (<3 epochs) were handled using cubic spline interpolation, totaling 45 points, while medium-term missing data (3-10 epochs) were handled using linear interpolation, totaling 15 points.

[0079] During the segmented discontinuity detection phase, two measurement arc segments were detected, with the arc segment interval being a natural interval (caused by the satellite entering and leaving the visible range of the ground station), and the interval was kept unchanged without interpolation. Anomaly segment detection employed a multi-method fusion strategy based on gradient detection, Pettitt test, and Mann-Kendall test, but no abnormal breakpoints were detected.

[0080] During the adaptive algorithm selection phase, a signal-to-noise ratio (SNR) of 42 dB (indicating high-quality data) was obtained through data feature analysis, along with data smoothing metrics. =0.008 (indicating very smooth data), and the non-stationarity detection sliding window variance change rate Δσ² / σ² = 0.08 < 0.5 (indicating stationarity). According to the decision tree rules, since SNR = 42 dB > 40 dB and the data is smooth, the Savitzky-Golay filtering algorithm is selected, with a recommended window size w = 15 and polynomial order n = 4. After grid search parameter optimization, the optimal parameters are determined to be window = 15 and poly = 4.

[0081] In the multidimensional quality assessment phase, the integrity score C = 100 × 4992 / 5000 = 99.84 points; in the accuracy assessment, RMSE = 0.08 m, Range = 300 m, then the accuracy score A = 100 × (1 0.08 / 300) = 99.97 points; the number of outliers in the consistency assessment is N_outlier = 8, then the consistency score Co = 100 × (1 The timeliness score is calculated as follows: 8 / 4992) = 99.84 points; In the timeliness assessment, the data delay Δt = 5s, τ = 3600s, then the timeliness score T = 100 × exp(-5 / 3600) = 99.86 points. Through weighted processing, the overall quality score is obtained. It scores 99.90 points, which is excellent (90-100 points), and will be output directly.

[0082] That is, after the original 5,000 data points went through a complete data preprocessing process, 4,992 high-quality data points were output, with a comprehensive data quality score of 99.90, which meets the requirements for high-precision orbit determination.

[0083] This embodiment employs the exact same preprocessing framework as the previous embodiment to process GNSS auxiliary data (RINEX 4.02 format, 0.1Hz sampling rate). It verifies the broad applicability of the method to different data sources and sampling rates, and also confirms its versatility and robustness under different data formats (RINEX 4.02) and different data characteristics (higher SNR 42dB). For high-quality data, it adaptively selects a larger window (w=15) and a higher polynomial order (n=4), fully demonstrating the algorithm's and parameters' adaptive capabilities and effectively ensuring the accuracy and reliability of GNSS auxiliary data preprocessing. Furthermore, the system automatically identifies the RINEX 4.02 format and calls the corresponding parser through a strategy pattern, without requiring modification of the underlying code, showcasing the advantages of a plug-in architecture.

[0084] As can be seen from the two embodiments above, the medium Earth orbit (MEO) satellite orbit measurement data preprocessing method described in this application achieves pipelined, orderly processing through hierarchical, progressive data cleaning; improves processing efficiency and detection accuracy through multi-dimensional parallel detection; overcomes the limitations of fixed algorithms through adaptive algorithm selection; preserves the physical meaning of data through intelligent identification of segmented discontinuities; and ensures high-quality output data through multi-dimensional quality assessment and graded processing. Compared with traditional preprocessing methods based on fixed processes, the method of this application achieves significant improvements in data quality pass rate, processing efficiency, and system scalability, and can meet the engineering requirements for near real-time, high-precision preprocessing of MEO satellite orbit measurement data.

[0085] It should be understood that although the steps in the flowcharts of the embodiments described above are shown sequentially according to the arrows, these steps are not necessarily executed in the order indicated by the arrows. Unless explicitly stated herein, there is no strict order restriction on the execution of these steps, and they can be executed in other orders. Moreover, at least some steps in the flowcharts of the embodiments described above may include multiple steps or multiple stages. These steps or stages are not necessarily completed at the same time, but can be executed at different times. The execution order of these steps or stages is not necessarily sequential, but can be performed alternately or in turn with other steps or at least some of the steps or stages in other steps. It is understood that the steps in different embodiments can be freely combined as needed, and all non-contradictory solutions formed by such combinations are within the scope of protection of this application.

[0086] Based on the same inventive concept, this application also provides an apparatus for implementing the aforementioned method for preprocessing medium Earth orbit measurement data. The solution provided by this apparatus is similar to the solution described in the above method. Therefore, the specific limitations in one or more embodiments of the medium Earth orbit measurement data preprocessing apparatus provided below can be found in the limitations of the medium Earth orbit measurement data preprocessing method described above, and will not be repeated here.

[0087] In one exemplary embodiment, such as Figure 3 As shown, a preprocessing device for medium Earth orbit satellite orbit measurement data is provided, comprising: a data access module 302, a data cleaning module 304, an arc segment identification module 306, a data filtering module 308, and a data output module 310, wherein: The data access module 302 is used to perform format recognition, integrity check and timestamp verification on the acquired raw orbit measurement data to obtain verified measurement data; the raw orbit measurement data includes at least one of X-band two-way ranging data, velocity data and GNSS auxiliary data from medium orbit relay satellites; The data cleaning module 304 is used to perform data cleaning processing on the verified measurement data to obtain cleaned measurement data; the data cleaning processing includes sequentially performing the measurement data in the order of time series sorting, duplicate data detection, outlier detection and missing value handling. The arc segment identification module 306 is used to identify the boundary of the measurement arc segment based on the time interval threshold, divide the cleaned measurement data into at least one measurement arc segment, and verify the division result of the measurement arc segment by combining the ground station position and orbit parameters to calculate the theoretical visibility, so as to distinguish the natural arc segment interval between adjacent measurement arc segments and the abnormal segment discontinuity within the measurement arc segment, and take corresponding processing strategies. The data filtering module 308 is used to perform multi-dimensional feature analysis of the processed measurement data in the time domain, frequency domain, and statistical domain to obtain data features, and to determine the filtering algorithm and resampling algorithm based on the data features, so as to filter, denoise, and resample the processed measurement data. The data output module 310 is used to calculate the quality score of the resampled measurement data from four dimensions: data integrity, accuracy, consistency and timeliness. The quality scores of each dimension are weighted and combined to obtain a comprehensive quality score. The data is then graded and processed according to the comprehensive quality score, and the preprocessed data is output.

[0088] In one exemplary embodiment, the data access module is specifically used for: The file header magic number of the original track measurement data is detected to identify the data format and version, and the parser corresponding to the data format and version is called to parse the original track measurement data; The parsed orbit measurement data is subjected to an integrity check, which includes three levels of checks: format level, semantic level, and reference level. The timestamps of the orbital measurement data that have passed the integrity check are unified to the unified time system of GPS or UTC. Leap second correction is performed based on the cumulative leap second difference between GPS and UTC, and the time interval between adjacent epochs is checked to ensure that it is consistent with the nominal sampling interval in order to obtain the verified measurement data.

[0089] In one exemplary embodiment, the duplicate data detection includes: The measurement data sorted according to time series is subjected to parallel duplicate detection from three dimensions: time, space, and attribute. The time dimension is used to detect whether the time difference between any two data records is less than a preset time threshold; the space dimension is used to detect whether the spatial distance between the two data records in a unified coordinate system is less than a preset distance threshold; and the attribute dimension is used to detect whether the difference between the pseudocode ranging observation values ​​of the two data records is less than a preset attribute threshold. The detection results from the three dimensions are aggregated through weighted voting to determine the duplicate data group, and the data record with the minimum uncertainty is retained within the duplicate data group.

[0090] In one exemplary embodiment, the outlier detection includes: Three statistical methods—Z-score, median absolute deviation, and interquartile range—were used to perform parallel anomaly detection on the observation sequence of the measurement data after duplicate data detection, in order to determine the anomaly detection result for each sequence point. Physical constraints are verified for each sequence point based on the physical constraints of medium Earth orbit satellite orbit dynamics to determine the anomaly verification result for each sequence point; the physical constraints include position constraints, velocity constraints, and acceleration constraints. The observation sequence is checked for temporal continuity to determine the temporal continuity result for each sequence point. The anomaly detection result, anomaly verification result, and temporal continuity result of each sequence point are weighted and voted on according to preset weights. When the weighted score exceeds the preset score threshold, the corresponding sequence point is determined to be an anomaly.

[0091] In one exemplary embodiment, the missing value processing includes: A theoretical epoch sequence is generated based on the start and end times of the measurement data and the nominal sampling interval. The theoretical epoch sequence is matched with the measured epochs one by one, and the unmatched theoretical epochs are identified as missing data. When the length of the missing interval is less than the first preset number of epochs, cubic spline interpolation is used to interpolate the missing interval to obtain the interpolated epoch sequence. When the length of the missing interval is not less than the first preset epoch number and not greater than the second preset epoch number, the missing interval is interpolated by linear interpolation or orbit extrapolation to obtain the interpolated epoch sequence. When the length of the missing interval is greater than the second preset epoch number, the missing interval is marked as an unusable segment and no interpolation is performed.

[0092] In an exemplary embodiment, the detection of abnormal segment discontinuities employs a multi-method fusion strategy that combines gradient detection, Pettitt test, and Mann-Kendall test. For identified discontinuous segments, a segmented independent processing, smooth transition, or label preservation strategy is adopted; for identified natural arc segment intervals, no interpolation is performed.

[0093] In one exemplary embodiment, the data filtering module is further configured to: The noise level is estimated based on the energy of the high-frequency components of the measurement data, and the signal-to-noise ratio is calculated; a normalized smoothness metric is calculated based on the second-order difference of the observation sequence; periodicity detection and non-stationarity detection are performed, and the non-stationarity detection determines whether the measurement data is stationary by calculating the rate of change of the variance of the sliding window. When the signal-to-noise ratio is greater than the preset signal-to-noise ratio threshold and the smoothness metric is less than the preset smoothness threshold, the filtering algorithm is determined to be Savitzky-Golay filtering. When the signal-to-noise ratio is not greater than the preset signal-to-noise ratio threshold, the filtering algorithm is determined to be moving average filtering; when the data has outlier characteristics, the filtering algorithm is determined to be median filtering. When the proportion of high-frequency components is greater than the preset proportion threshold, the filtering algorithm is determined to be low-pass filtering.

[0094] In an exemplary embodiment, the weighted combination of the quality scores of each dimension to obtain the comprehensive quality score includes: Where C represents the completeness score and A represents the accuracy score. For consistency scoring, T represents timeliness scoring. , , and For the corresponding weights, This is the overall quality score.

[0095] Each module in the aforementioned medium-Earth orbit satellite orbit measurement data preprocessing device can be implemented entirely or partially through software, hardware, or a combination thereof. These modules can be embedded in or independent of the processor in a computer device, or stored in the computer device's memory as software, so that the processor can call and execute the corresponding operations of each module.

[0096] In one exemplary embodiment, a computer device is provided, the internal structure of which can be as shown in the figure. Figure 4 As shown, the computer device includes a processor, memory, input / output interface, communication interface, display unit, and input device. The processor, memory, and input / output interface are connected via a system bus, and the communication interface, display unit, and input device are also connected to the system bus via the input / output interface. The processor provides computational and control capabilities. The memory includes non-volatile storage media and internal memory. The non-volatile storage media stores the operating system and computer programs. The internal memory provides an environment for the operation of the operating system and computer programs stored in the non-volatile storage media. The input / output interface is used for exchanging information between the processor and external devices. The communication interface is used for wired or wireless communication with external terminals; wireless communication can be achieved through Wi-Fi, mobile cellular networks, Near Field Communication (NFC), or other technologies. When executed by the processor, the computer program implements a method for preprocessing medium-Earth orbit satellite orbit measurement data. The display unit is used to form a visually visible image and can be a display screen, projection device, or virtual reality imaging device. The display screen can be an LCD screen or an e-ink screen. The input device of the computer device can be a touch layer covering the display screen, or buttons, trackballs, or touchpads set on the casing of the computer device, or external keyboards, touchpads, or mice, etc.

[0097] Those skilled in the art will understand that Figure 4 The structure shown is merely a block diagram of a portion of the structure related to the present application and does not constitute a limitation on the computer device to which the present application is applied. Specific computer devices may include more or fewer components than those shown in the figure, or combine certain components, or have different component arrangements.

[0098] In one exemplary embodiment, a computer device is provided, including a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program to implement the steps in the above-described method embodiments.

[0099] In one embodiment, a computer-readable storage medium is provided having a computer program stored thereon, which, when executed by a processor, implements the steps in the above method embodiments.

[0100] In one embodiment, a computer program product is provided, including a computer program that, when executed by a processor, implements the steps in the above method embodiments.

[0101] It should be noted that the user information (including but not limited to user device information, user personal information, etc.) and data (including but not limited to data used for analysis, data stored, data displayed, etc.) involved in this application are all information and data authorized by the user or fully authorized by all parties, and the collection, use and processing of the relevant data must comply with relevant regulations.

[0102] Those skilled in the art will understand that all or part of the processes in the methods of the above embodiments can be implemented by a computer program instructing related hardware. The computer program can be stored in a non-volatile computer-readable storage medium, and when executed, it can include the processes of the embodiments of the above methods. Any references to memory, databases, or other media used in the embodiments provided in this application can include at least one of non-volatile memory and volatile memory. Non-volatile memory can include read-only memory (ROM), magnetic tape, floppy disk, flash memory, optical memory, high-density embedded non-volatile memory, resistive random access memory (ReRAM), magnetic random access memory (MRAM), ferroelectric random access memory (FRAM), phase change memory (PCM), graphene memory, etc. Volatile memory can include random access memory (RAM) or external cache memory, etc. By way of illustration and not limitation, RAM can take many forms, such as Static Random Access Memory (SRAM) or Dynamic Random Access Memory (DRAM). The databases involved in the embodiments provided in this application may include at least one type of relational database and non-relational database. Non-relational databases may include, but are not limited to, blockchain-based distributed databases. The processors involved in the embodiments provided in this application may be general-purpose processors, central processing units, graphics processing units, digital signal processors, programmable logic devices, quantum computing-based data processing logic devices, artificial intelligence (AI) processors, etc., and are not limited to these.

[0103] The technical features of the above embodiments can be combined in any way. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this application.

[0104] The embodiments described above are merely illustrative of several implementation methods of this application, and while the descriptions are specific and detailed, they should not be construed as limiting the scope of this patent application. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of this application, and these all fall within the protection scope of this application. Therefore, the protection scope of this application should be determined by the appended claims.

Claims

1. A method for preprocessing orbit measurement data of medium Earth orbit satellites, characterized in that, The method includes: The acquired raw orbit measurement data is subjected to format recognition, integrity check and timestamp verification to obtain verified measurement data; the raw orbit measurement data includes at least one of X-band two-way ranging data, velocity data and GNSS auxiliary data from medium Earth orbit relay satellites; The validated measurement data is cleaned to obtain cleaned measurement data. The data cleaning process includes sequentially sorting the measurement data according to time series, detecting duplicate data, detecting outliers, and handling missing values. The boundary of the measurement arc segment is identified based on the time interval threshold. The cleaned measurement data is divided into at least one measurement arc segment. The theoretical visibility of the measurement arc segment is calculated in combination with the ground station location and orbit parameters to verify the division result of the measurement arc segment. This is to distinguish between the natural arc segment interval between adjacent measurement arc segments and the abnormal segment discontinuity within the measurement arc segment, and to take corresponding processing strategies. Multidimensional feature analysis in the time domain, frequency domain, and statistical domain is performed on the processed measurement data to obtain data features. Based on the data features, filtering algorithms and resampling algorithms are determined to filter, denoise, and resample the processed measurement data. The quality scores of the resampled measurement data are calculated from four dimensions: completeness, accuracy, consistency, and timeliness. The quality scores of each dimension are weighted and combined to obtain a comprehensive quality score. The data is then graded and processed according to the comprehensive quality score, and the preprocessed data is output.

2. The method according to claim 1, characterized in that, The process of performing format recognition, integrity checks, and timestamp verification on the acquired raw orbit measurement data to obtain verified measurement data includes: The file header magic number of the original track measurement data is detected to identify the data format and version, and the parser corresponding to the data format and version is called to parse the original track measurement data; The parsed orbit measurement data is subjected to an integrity check, which includes three levels of checks: format level, semantic level, and reference level. The timestamps of the orbital measurement data that have passed the integrity check are unified to the unified time system of GPS or UTC. Leap second correction is performed based on the cumulative leap second difference between GPS and UTC, and the time interval between adjacent epochs is checked to ensure that it is consistent with the nominal sampling interval in order to obtain the verified measurement data.

3. The method according to claim 1, characterized in that, The duplicate data detection includes: The measurement data sorted according to time series is subjected to parallel duplicate detection from three dimensions: time, space, and attribute. The time dimension is used to detect whether the time difference between any two data records is less than a preset time threshold; the space dimension is used to detect whether the spatial distance between the two data records in a unified coordinate system is less than a preset distance threshold; and the attribute dimension is used to detect whether the difference between the pseudocode ranging observation values ​​of the two data records is less than a preset attribute threshold. The detection results from the three dimensions are aggregated through weighted voting to determine the duplicate data group, and the data record with the minimum uncertainty is retained within the duplicate data group.

4. The method according to claim 3, characterized in that, The outlier detection includes: Three statistical methods—Z-score, median absolute deviation, and interquartile range—were used to perform parallel anomaly detection on the observation sequence of the measurement data after duplicate data detection, in order to determine the anomaly detection result for each sequence point. Physical constraints are verified for each sequence point based on the physical constraints of medium Earth orbit satellite orbit dynamics to determine the anomaly verification result for each sequence point; the physical constraints include position constraints, velocity constraints, and acceleration constraints. The observation sequence is checked for temporal continuity to determine the temporal continuity result for each sequence point. The anomaly detection result, anomaly verification result, and temporal continuity result of each sequence point are weighted and voted on according to preset weights. When the weighted score exceeds the preset score threshold, the corresponding sequence point is determined to be an anomaly.

5. The method according to claim 4, characterized in that, The missing value processing includes: A theoretical epoch sequence is generated based on the start and end times of the measurement data and the nominal sampling interval. The theoretical epoch sequence is matched with the measured epochs one by one, and the unmatched theoretical epochs are identified as missing data. When the length of the missing interval is less than the first preset number of epochs, cubic spline interpolation is used to interpolate the missing interval to obtain the interpolated epoch sequence. When the length of the missing interval is not less than the first preset epoch number and not greater than the second preset epoch number, the missing interval is interpolated by linear interpolation or orbit extrapolation to obtain the interpolated epoch sequence. When the length of the missing interval is greater than the second preset epoch number, the missing interval is marked as an unusable segment and no interpolation is performed.

6. The method according to claim 1, characterized in that, The detection of abnormal segment discontinuities adopts a multi-method fusion strategy that combines gradient detection, Pettitt test and Mann-Kendall test; For identified discontinuous segments, a segmented independent processing, smooth transition, or label preservation strategy is adopted; for identified natural arc segment intervals, no interpolation is performed.

7. The method according to any one of claims 1 to 6, characterized in that, The process involves performing multidimensional feature analysis on the processed measurement data in the time domain, frequency domain, and statistical domain to obtain data features, and determining a filtering algorithm based on these data features, including: The noise level is estimated based on the energy of the high-frequency components of the measurement data, and the signal-to-noise ratio is calculated; a normalized smoothness metric is calculated based on the second-order difference of the observation sequence; periodicity detection and non-stationarity detection are performed, and the non-stationarity detection determines whether the measurement data is stationary by calculating the rate of change of the variance of the sliding window. When the signal-to-noise ratio is greater than the preset signal-to-noise ratio threshold and the smoothness metric is less than the preset smoothness threshold, the filtering algorithm is determined to be Savitzky-Golay filtering. When the signal-to-noise ratio is not greater than the preset signal-to-noise ratio threshold, the filtering algorithm is determined to be moving average filtering; when the data has outlier characteristics, the filtering algorithm is determined to be median filtering. When the proportion of high-frequency components is greater than the preset proportion threshold, the filtering algorithm is determined to be low-pass filtering.

8. The method according to any one of claims 1 to 6, characterized in that, The weighted combination of quality scores from each dimension to obtain the overall quality score includes: Where C represents the completeness score and A represents the accuracy score. For consistency scoring, T represents timeliness scoring. , , and For the corresponding weights, This is the overall quality score.

9. A preprocessing device for medium Earth orbit satellite orbit measurement data, characterized in that, The device includes: The data access module is used to perform format recognition, integrity check, and timestamp verification on the acquired raw orbit measurement data to obtain verified measurement data; the raw orbit measurement data includes at least one of X-band two-way ranging data, velocity data, and GNSS auxiliary data from medium Earth orbit relay satellites; The data cleaning module is used to clean the validated measurement data to obtain cleaned measurement data. The data cleaning process includes sequentially performing the measurement data sorting according to time series, duplicate data detection, outlier detection, and missing value handling. The arc segment identification module is used to identify the boundary of the measurement arc segment based on the time interval threshold, divide the cleaned measurement data into at least one measurement arc segment, and verify the division result of the measurement arc segment by combining the ground station position and orbit parameters to calculate the theoretical visibility, so as to distinguish the natural arc segment interval between adjacent measurement arc segments and the abnormal segment discontinuity within the measurement arc segment, and take corresponding processing strategies. The data filtering module is used to perform multi-dimensional feature analysis of the processed measurement data in the time domain, frequency domain, and statistical domain to obtain data features. Based on the data features, the filtering algorithm and resampling algorithm are determined to filter, denoise, and resample the processed measurement data. The data output module is used to calculate the quality score of the resampled measurement data from four dimensions: completeness, accuracy, consistency and timeliness. The quality scores of each dimension are weighted and combined to obtain the comprehensive quality score. The data is then graded and processed according to the comprehensive quality score, and the preprocessed data is output.

10. A computer device comprising a memory and a processor, wherein the memory stores a computer program, characterized in that, When the processor executes the computer program, it implements the steps of the method according to any one of claims 1 to 8.

11. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by a processor, it implements the steps of the method according to any one of claims 1 to 8.

12. A computer program product, comprising a computer program, characterized in that, When the computer program is executed by a processor, it implements the steps of the method according to any one of claims 1 to 8.

Images