Multi-source abnormal information real-time detection method for real-time orbit determination of Beidou system

Abstract

The invention discloses a multi-source abnormal information real-time detection method for real-time orbit determination of a Beidou system, and relates to the field of satellite abnormal detection. According to the method, information obtained through real-time orbit determination filtering calculation is used as a detection factor, classification selection of the detection factor is carried out through a phase jump detection threshold value, and the real-time orbit determination of the Beidou system is realized; according to the method, the anomaly detection factors are cooperatively screened through the isolation forest algorithm and the DBSCAN clustering analysis method, the anomaly generation source is judged through the proportion of the anomaly detection factors, anomaly of different magnitudes, different characteristics and different sources can be effectively detected, and the influence of anomaly information on the service capacity of the Beidou system is greatly reduced.

Description

Technical Field

[0001] The present invention relates to the field of satellite anomaly detection, and particularly to a real-time detection method for multi-source anomaly information in the real-time orbit determination of the Beidou system. Background Art

[0002] During the on-orbit operation of Beidou satellites, they are faced with an anomaly information flow that is rich in sources, diverse in characteristics, and dynamically changing over time. Its existence poses a serious threat to the high-precision service capabilities of the Beidou system. Therefore, clarifying the specific sources and manifestation forms of multi-source anomaly information, constructing detection factors that can sensitively respond to multi-source anomaly situations, and studying the modeling algorithms of detection factors are of great significance for the stable operation of the Beidou system and maintaining real-time high-precision services.

[0003] When a hardware failure occurs in the signal transponder of a Beidou satellite, it will cause an anomaly in the transmitted signal, thereby generating gross error observations (observation anomalies); when GEO satellites and IGSO satellites are in long-term on-orbit operation, they are affected by disturbing forces such as the non-spherical gravitational force of the earth, and will gradually deviate from the preset orbit. In order to maintain the orbit position, orbit maneuvers will occur frequently (state anomalies); affected by factors such as the atomic clock manufacturing principle, active phase modulation and frequency modulation by the master control station, and changes in temperature and humidity, the atomic clocks of Beidou satellites will experience phase jumps and frequency jumps (state anomalies). The unknown, sudden, and highly outlier observation anomalies and state anomalies pose a major challenge to the reliability and integrity of the Beidou system. Summary of the Invention

[0004] Aiming at the above deficiencies in the prior art, the real-time detection method for multi-source anomaly information in the real-time orbit determination of the Beidou system provided by the present invention can quickly detect the observation anomalies and state anomalies of Beidou satellites.

[0005] In order to achieve the above invention purpose, the technical solution adopted by the present invention is as follows:

[0006] Provide a real-time detection method for multi-source anomaly information in the real-time orbit determination of the Beidou system, which includes the following steps:

[0007] S1. Calculate the predicted state value of the satellite within the current epoch by obtaining the state transition matrix between the current epoch and the previous epoch of the satellite;

[0008] S2. Obtain the filtering innovation value based on the predicted state value of the satellite and the observed value of the satellite;

[0009] S3. Set the phase jump detection threshold, and determine whether there is a filtering innovation value greater than or equal to the phase jump detection threshold within the epoch to be detected. If so, use the filtering innovation value greater than or equal to the phase jump detection threshold as the detection factor; otherwise, use all the filtering innovation values within the epoch to be detected as the detection factor;

[0010] S4. Calculate the outlier score corresponding to the detection factor through the Isolation Forest algorithm, and mark the detection factors with outlier scores greater than or equal to the score threshold as abnormal detection factors;

[0011] S5. Cluster the detection factors through the DBSCAN clustering analysis method to obtain abnormal detection factors;

[0012] S6. Determine the source of the anomaly by calculating the proportion of the abnormal detection factors, and complete the real-time detection of multi-source anomaly information for the real-time orbit determination of the Beidou system.

[0013] Further, the specific method for obtaining the state transition matrix of the satellite between the current epoch and the previous epoch in step S1 is as follows:

[0014] Based on the earth perturbation force, solar and lunar gravitational perturbation, earth tide perturbation, relativistic effect, and solar radiation pressure perturbation acting on the GEO satellites, IGSO satellites, and MEO satellites in the Beidou system, establish the perturbation motion equation, and use numerical methods to solve for the state transition matrix Φ of the satellite between the current epoch and the previous epoch k-1,k ;

[0015] The expression for the predicted state value of the satellite within the current epoch in step S1 is:

[0016]

[0017] where is the predicted state value of the satellite within the current epoch; X k-1 is the satellite state parameter within the previous epoch.

[0018] Further, the calculation expression for the filter innovation value is:

[0019]

[0020] where is the filter innovation value matrix within the current epoch, is the m-th filter innovation value in the current epoch; A k is the design matrix; L k is the observation value matrix of the satellite within the current epoch.

[0021] Further, the specific method for calculating the outlier corresponding to the detection factor through the Isolation Forest algorithm in step S4 includes the following sub-steps:

[0022] S4-1. Obtain a number of filter innovation values using the same method as in step S2. Starting from the epoch when the filter innovation values converge, select the filter innovation values within N epochs to form a training set, and train the Isolation Forest algorithm with the training set to obtain the trained Isolation Forest algorithm;

[0023] S4-2. Traverse each isolated tree of the trained isolation forest algorithm with the detection factor in the order of epoch sequence, and calculate the average height of the detection factor in the isolation forest;

[0024] S4-3. Normalize the average heights of all detection factors in the isolation forest, and calculate the outlier scores corresponding to each detection factor.

[0025] Furthermore, the specific method for training the isolation forest algorithm in step S4-1 includes the following sub-steps:

[0026] S4-1-1. Randomly select ψ filtered innovation values from the training set as a sample subset;

[0027] S4-1-2. For any sample subset, perform a splitting operation:

[0028] Randomly select a feature, and randomly select a splitting point within the value range of this feature to divide the current data into two parts; among them, the data less than the splitting point in the randomly selected feature is placed in the left child node of the current node, and the data greater than or equal to the splitting point is placed in the right child node of the current node; the data at the first splitting is the sample subset;

[0029] S4-1-3. Repeat the splitting operation in step S4-1-2 until the data in the sample subset cannot be further divided, obtaining leaf nodes, and thus completing the construction of an isolated tree;

[0030] S4-1-4. Repeat steps S4-1-1 to S4-1-3 until T isolated trees are generated, and thus obtain the trained isolation forest algorithm.

[0031] Furthermore, the value of ψ is 256, and the value of T is 100.

[0032] Furthermore, the expression for calculating the outlier score corresponding to each detection factor in step S4-3 is:

[0033] (ψ - 1) = ln(ψ - 1) + γ

[0034]

[0035]

[0036] where H(ψ - 1) is the harmonic number; ln(.) is the natural logarithm; γ is the Euler constant; c(ψ) is the average height of the nodes in the isolated tree; S(x, ψ) is the outlier score corresponding to the detection factor x; E(h(x)) is the height of the detection factor x after normalization; the score threshold in step S4 is 0.6.

[0037] Further, in step S5, the detection factors are clustered by the DBSCAN clustering analysis method, and the specific method for obtaining the abnormal detection factors includes the following sub-steps:

[0038] S4-1. Set the neighborhood radius R and the minimum number of points minpoints within the neighborhood;

[0039] S4-2. Traverse all the detection factors of the epoch to be detected. For the detection factors with the number of points within the radius greater than or equal to minpoints, include them as core points in the core point list; form a temporary clustering cluster corresponding to the core point for the detection factors within the radius of the core point;

[0040] S4-3. For each temporary clustering cluster, check whether the detection factors in it are core points. If so, merge the temporary clustering cluster corresponding to the detection factor as the core point with the current temporary clustering cluster until each detection factor in the temporary clustering cluster is not in the core point list, or all the detection factors within the neighborhood radius R of each detection factor in the temporary clustering cluster are in the same temporary clustering cluster, and obtain the clustering cluster;

[0041] S4-4. Regard the detection factors not classified into the clustering cluster as abnormal detection factors.

[0042] Further, in step S6, the specific method for determining the source of the anomaly by calculating the proportion of abnormal detection factors includes the following sub-steps:

[0043] S5-1. Calculate the proportion of abnormal detection factors, and its expression is:

[0044]

[0045] where is the proportion of abnormal detection factors; represents the total number of detection factors; represents the total number of abnormal detection factors;

[0046] S5-2. When there is a filtering innovation value greater than or equal to the phase jump detection threshold in the epoch to be detected, and the proportion of abnormal detection factors is 100%, it is determined that the corresponding satellite has a phase jump fault;

[0047] When there is no filtering innovation value greater than or equal to the phase jump detection threshold in the epoch to be detected, and the proportion of abnormal detection factors of each satellite in the Beidou system is less than 90%, it is determined as an observation anomaly;

[0048] When there is no filtering innovation value greater than or equal to the phase jump detection threshold in the epoch to be detected, and the proportion of abnormal detection factors of the satellite is greater than or equal to 90% and less than 100%, it is determined that the satellite is in an abnormal state, and it is determined that the anomaly is caused by the frequency jump of the satellite atomic clock;

[0049] When there is no filtering innovation value greater than or equal to the phase jump detection threshold in the epoch to be detected, and the proportion of the satellite's anomaly detection factor is 100%, it is determined that the satellite is in an abnormal state, and it is determined that the anomaly is caused by satellite orbit maneuvering.

[0050] The beneficial effects of the present invention are as follows: This method uses the innovation obtained from real-time orbit determination filtering calculation as the detection factor, classifies and selects the detection factor through the phase jump detection threshold, collaboratively screens the anomaly detection factor through the Isolation Forest algorithm and the DBSCAN clustering analysis method, and determines the source of the anomaly through the proportion of the anomaly detection factor, which can effectively detect anomalies of different magnitudes, different characteristics and different sources, and greatly reduce the impact of anomaly information on the service ability of the Beidou system. Description of the Drawings

[0051] Figure 1 It is a schematic flow diagram of this method. Detailed Embodiments

[0052] The following describes the detailed embodiments of the present invention to facilitate those skilled in the art of the present technology to understand the present invention. However, it should be clear that the present invention is not limited to the scope of the detailed embodiments. For those of ordinary skill in the art of the present technology, as long as various changes are within the spirit and scope of the present invention defined and determined by the appended claims, these changes are obvious, and all inventions and creations using the concept of the present invention are within the scope of protection.

[0053] As Figure 1 shown, the multi-source anomaly information real-time detection method for real-time orbit determination of the Beidou system includes the following steps:

[0054] S1. Calculate the predicted state value of the satellite in the current epoch by obtaining the state transition matrix between the current epoch and the previous epoch of the satellite;

[0055] S2. Obtain the filtering innovation value based on the predicted state value of the satellite and the observed value of the satellite;

[0056] S3. Set the phase jump detection threshold, and determine whether there is a filtering innovation value greater than or equal to the phase jump detection threshold in the epoch to be detected. If so, use the filtering innovation value greater than or equal to the phase jump detection threshold as the detection factor; otherwise, use all the filtering innovation values in the epoch to be detected as the detection factor; the current epoch is the epoch to be detected;

[0057] S4. Calculate the anomaly score corresponding to the detection factor through the Isolation Forest algorithm, and record the detection factor with an anomaly score greater than or equal to the score threshold as the anomaly detection factor;

[0058] S5. Cluster the detection factors through the DBSCAN clustering analysis method to obtain the anomaly detection factors;

[0059] S6. Determine the source of anomalies by calculating the proportion of anomaly detection factors, and complete the real-time detection of multi-source anomaly information for real-time orbit determination of the Beidou system.

[0060] The specific method for obtaining the state transition matrix of the satellite between the current epoch and the previous epoch in step S1 is as follows:

[0061] Based on the earth disturbing force, solar and lunar gravitational perturbations, earth tide perturbations, relativistic effects, and solar radiation pressure perturbations acting on the GEO satellites, IGSO satellites, and MEO satellites in the Beidou system, establish the perturbation motion equation, and use numerical methods to solve for the state transition matrix Φ between the current epoch and the previous epoch of the satellite. k-1,k ;

[0062] The expression for the predicted state value of the satellite within the current epoch in step S1 is:

[0063]

[0064] where is the predicted state value of the satellite within the current epoch; X k-1 is the satellite state parameter within the previous epoch.

[0065] The calculation expression for the filter innovation value is:

[0066]

[0067] where is the filter innovation value matrix within the current epoch, is the m-th filter innovation value in the current epoch; A k is the design matrix; L k is the observation value matrix of the satellite within the current epoch.

[0068] The specific method for calculating the anomaly value corresponding to the detection factor by the isolation forest algorithm in step S4 includes the following sub-steps:

[0069] S4-1. Obtain a number of filter innovation values using the same method as in step S2. Starting from the epoch when the filter innovation values converge, select the filter innovation values within N epochs to form a training set, and train the isolation forest algorithm using the training set to obtain the trained isolation forest algorithm;

[0070] S4-2. Let the detection factor traverse each isolated tree of the trained isolation forest algorithm in the order of epochs, and calculate the average height of the detection factor in the isolation forest;

[0071] S4-3. Normalize the average height of all detection factors in the isolation forest and calculate the outlier score corresponding to each detection factor.

[0072] The specific method for training the isolation forest algorithm in step S4-1 includes the following sub-steps:

[0073] S4-1-1. Randomly select ψ filtering innovation values from the training set as a sample subset;

[0074] S4-1-2. For any sample subset, perform a splitting operation:

[0075] Randomly select a feature and randomly select a splitting point within the value range of this feature to divide the current data into two parts; among them, the data less than the splitting point in the randomly selected feature is placed in the left child node of the current node, and the data greater than or equal to the splitting point is placed in the right child node of the current node; the data at the first splitting is the sample subset;

[0076] S4-1-3. Repeat the splitting operation in step S4-1-2 until the data in the sample subset cannot be further divided, obtain leaf nodes, and thus complete the construction of an isolation tree;

[0077] S4-1-4. Repeat steps S4-1-1 to S4-1-3 until T isolation trees are generated, and thus obtain the trained isolation forest algorithm.

[0078] The value of ψ is 256, and the value of T is 100.

[0079] The expression for calculating the outlier score corresponding to each detection factor in step S4-3 is:

[0080] (ψ - 1) = ln(ψ - 1) + γ

[0081]

[0082]

[0083] Where H(ψ - 1) is the harmonic number; ln(.) is the natural logarithm; γ is the Euler constant; c(ψ) is the average height of the nodes in the isolation tree; S(x, ψ) is the outlier score corresponding to the detection factor x; E(h(x)) is the height of the detection factor x after normalization; the score threshold in step S4 is 0.6.

[0084] 8. The real-time multi-source anomaly information real-time detection method for Beidou system real-time orbit determination according to claim 1, characterized in that in step S5, the detection factors are clustered by the DBSCAN clustering analysis method, and the specific method for obtaining the anomaly detection factors includes the following sub-steps:

[0085] S4-1. Set the neighborhood radius R and the minimum number of points minpoints in the neighborhood;

[0086] S4-2. Traverse all the detection factors of the epoch to be detected. For the detection factors whose number of points within the radius is greater than or equal to minpoints, include them as core points in the core point list; form a temporary clustering cluster corresponding to the core point for the detection factors within the radius of the core point;

[0087] S4-3. For each temporary clustering cluster, check whether the detection factors in it are core points. If so, merge the temporary clustering cluster corresponding to the detection factor as the core point with the current temporary clustering cluster until each detection factor in the temporary clustering cluster is not in the core point list, or all the detection factors within the neighborhood radius R of each detection factor in the temporary clustering cluster are in the same temporary clustering cluster, to obtain the clustering cluster;

[0088] S4-4. Regard the detection factors not classified into the clustering cluster as abnormal detection factors.

[0089] In this embodiment, the values of both R and minpoints are 4.

[0090] The specific method for determining the source of the anomaly by calculating the proportion of abnormal detection factors in step S6 includes the following sub-steps:

[0091] S5-1. Calculate the proportion of abnormal detection factors, and its expression is:

[0092]

[0093] where is the proportion of abnormal detection factors; represents the total number of detection factors; represents the total number of abnormal detection factors;

[0094] S5-2. When there is a filtering innovation value greater than or equal to the phase jump detection threshold in the epoch to be detected, and the proportion of abnormal detection factors is 100%, it is determined that the corresponding satellite has a phase jump fault;

[0095] When there is no filtering innovation value greater than or equal to the phase jump detection threshold in the epoch to be detected, and the proportion of abnormal detection factors of each satellite in the Beidou system is less than 90%, it is determined as an observation anomaly;

[0096] When there is no filtering innovation value greater than or equal to the phase jump detection threshold in the epoch to be detected, and the proportion of abnormal detection factors of the satellite is greater than or equal to 90% and less than 100%, it is determined that the satellite is in an abnormal state, and it is determined that the anomaly is caused by the frequency jump of the satellite atomic clock;

[0097] When there is no filtering innovation value greater than or equal to the phase jump detection threshold within the epoch to be detected, and the proportion of the satellite's anomaly detection factor is 100%, it is determined that the satellite is in an abnormal state, and it is determined that the anomaly is caused by satellite orbit maneuvering.

[0098] In the specific implementation process, taking the anomaly detection result of epoch 102 as an example, k is the satellite PRN number, i is the number of abnormal innovations of the satellite within this epoch; m is the total number of innovations of the satellite within this epoch, and P is the abnormal proportion. The detection results by the Isolation Forest algorithm are shown in Table 1 below.

[0099] Table 1

[0100]

[0101]

[0102] The detection results by the DBSCAN clustering analysis method are shown in Table 2 below.

[0103] Table 2

[0104]

[0105]

[0106] During the anomaly detection of epoch 102 observations, due to the influence of the visibility conditions of the simulated satellites or the failure to assign gross errors to some satellites when generating simulated gross errors, some satellites have no detection results.

[0107] In the specific implementation process, for the anomaly detection of a certain epoch, k is the satellite PRN number, i is the number of abnormal innovations of the satellite within this epoch; m is the total number of innovations of the satellite within this epoch, and P is the abnormal proportion. The detection results by the Isolation Forest algorithm are shown in Table 3 below.

[0108] Table 3

[0109]

[0110]

[0111]

[0112] The detection results by the DBSCAN clustering analysis method are shown in Table 4 below.

[0113]

[0114]

[0115] As can be seen from Table 3 and Table 4, PRN1 is the abnormal satellite that has maneuvered.

[0116] In the specific implementation process, the frequency jump data of the simulated satellite atomic clock is shown in Table 5 below.

[0117] Table 5

[0118]

[0119] The results of detecting PRN1&23 using the Isolation Forest algorithm are shown in Table 6 below.

[0120] Table 6

[0121]

[0122]

[0123] The results of detecting PRN7&22 using the Isolation Forest algorithm are shown in Table 7 below.

[0124] Table 7

[0125]

[0126]

[0127] The results of detecting PRN1&23 using the DBSCAN clustering analysis method are shown in Table 8 below.

[0128] Table 8

[0129] k i m P 1 11 33 0.333 2 1 34 0.029 3 1 35 0.029 4 1 34 0.029 5 1 33 0.030 6 1 33 0.030 7 2 40 0.050 8 2 43 0.047 9 1 41 0.024 11 1 40 0.025 12 1 41 0.024 14 2 44 0.045 15 1 44 0.023 16 1 42 0.024 17 1 42 0.024 18 1 41 0.024 19 1 42 0.024 20 2 40 0.050 21 1 37 0.027 22 1 45 0.022 23 22 39 0.564 24 1 42 0.024 25 1 38 0.026 26 1 42 0.024 27 1 40 0.025 28 1 38 0.026 29 2 42 0.048 30 1 42 0.024

[0130] The results of detecting PRN7&22 using the DBSCAN clustering analysis method are shown in Table 9 below.

[0131] Table 9

[0132]

[0133]

[0134] The statistical results of the frequency jump of PRN1 are shown in Table 10 below.

[0135] Table 10

[0136] Evaluation DBSCAN Isolation Forest Accuracy 0.9249 0.9590 Precision 1 0.7561 Recall 0.3333 0.9394 F1 Score 0.5000 0.8378

[0137] The statistical results of the frequency jump of PRN7 are shown in Table 11 below.

[0138] Table 11

[0139]

[0140]

[0141] As can be seen from the above experimental data, the present method uses the innovation obtained from real-time orbit determination filtering solution as the detection factor, classifies and selects the detection factor through the phase jump detection threshold. Both the DBSCAN clustering analysis method and the isolation forest algorithm can effectively screen out anomalies in the detection factor. By determining the proportion of the results of the collaborative screening of the two, the source of anomaly generation can be determined, and anomalies of different magnitudes, different characteristics and different sources can be effectively detected, greatly reducing the impact of anomaly information on the service ability of the Beidou system.

Claims

1. A method for real-time detection of multi-source anomaly information in real-time orbit determination of Beidou system, characterized in that: The following steps are involved: S1, calculating the predicted state value of the satellite in the current epoch by obtaining the state transfer matrix of the satellite between the current epoch and the previous epoch; S2. Obtaining a filtering innovation value based on the satellite predicted state value and the satellite observation value; S3, setting a phase jump detection threshold, determining whether there is a filtering innovation value greater than or equal to the phase jump detection threshold in the epoch to be detected, and if so, using the filtering innovation value greater than or equal to the phase jump detection threshold as a detection factor; Otherwise, all filtered innovation values ​​of the epoch to be detected are used as detection factors; S4. Calculate the outlier score corresponding to the detection factor by using the isolation forest algorithm, and record the detection factor whose outlier score is greater than or equal to the score threshold as the outlier detection factor; S5. Cluster the detection factors using the DBSCAN cluster analysis method to obtain anomaly detection factors; S6. Determine the source of the anomaly by calculating the proportion of anomaly detection factors, and complete the real-time detection of multi-source anomaly information in the real-time orbit determination of the Beidou system.

2. The method for real-time detection of multi-source abnormal information for Beidou system real-time orbit determination according to claim 1 is characterized in that: The specific method for obtaining the state transfer matrix of the satellite between the current epoch and the previous epoch in step S1 is: According to the earth perturbation, solar-moon gravitational perturbation, earth tidal perturbation, relativistic effect and solar light pressure perturbation on the GEO, IGSO and MEO satellites of the BeiDou system, the perturbation motion equation is established, and the state transfer matrix Φ between the current epoch and the previous epoch is obtained by numerical method. k-1,k ; The expression of the satellite prediction state value in the current epoch in step S1 is: in is the predicted satellite status value in the current epoch; X k-1 is the satellite status parameter in the previous epoch.

3. The method for real-time detection of multi-source abnormal information for Beidou system real-time orbit determination according to claim 2 is characterized in that: The calculation expression of the filter innovation value is: in is the filter innovation value matrix in the current epoch, is the mth filtering innovation value in the current epoch; A k is the design matrix; L k is the satellite observation matrix in the current epoch.

4. The method for real-time detection of multi-source abnormal information for Beidou system real-time orbit determination according to claim 1 is characterized in that: The specific method of calculating the outlier value corresponding to the detection factor by using the isolation forest algorithm in step S4 includes the following sub-steps: S4-1, using the same method as step S2 to obtain a number of filter innovation values, taking the epoch when the filter innovation value converges as the starting point, selecting the filter innovation values ​​within N epochs to form a training set, and training the isolation forest algorithm with the training set to obtain a trained isolation forest algorithm; S4-2, traversing each isolated tree of the trained isolation forest algorithm by the detection factor in an epoch-by-epoch sequence, and calculating the average height of the detection factor in the isolation forest; S4-3. Normalize the average height of all detection factors in the isolation forest and calculate the outlier score corresponding to each detection factor.

5. The method for real-time detection of multi-source abnormal information in Beidou system real-time orbit determination according to claim 4 is characterized in that: The specific method for training the isolation forest algorithm in step S4-1 includes the following sub-steps: S4-1-1, randomly select ψ filtered innovation values ​​from the training set as a sample subset; S4-1-2 For any sample subset, perform segmentation operation: Randomly select a feature, and randomly select a split point within the value range of the feature to divide the current data into two parts; the data less than the split point in the randomly selected feature is placed in the left child node of the current node, and the data greater than or equal to the split point is placed in the right child node of the current node; the data at the first split is the sample subset; S4-1-3, repeat the splitting operation of step S4-1-2 until the data of the sample subset cannot be divided any further, and a leaf node is obtained, thereby completing the construction of an isolated tree; S4-1-4, repeat steps S4-1-1 to S4-1-3 until T isolated trees are generated, thereby obtaining the trained isolation forest algorithm.

6. The method for real-time detection of multi-source abnormal information in Beidou system real-time orbit determination according to claim 5 is characterized in that: The value of ψ is 256 and the value of T is 100.

7. The method for real-time detection of multi-source abnormal information in Beidou system real-time orbit determination according to claim 5, characterized in that: The expression for calculating the outlier score corresponding to each detection factor in step S4-3 is: H(ψ-1)=ln(ψ-1)+γ Where H(ψ-1) is the harmonic number; ln(.) is the natural logarithm; γ is the Euler constant; c(ψ) is the average height of the nodes in the isolated tree; S(x,ψ) is the outlier score corresponding to the detection factor x; E(h(x)) is the normalized height of the detection factor x; and the score threshold in step S4 is 0.

6.

8. The method for real-time detection of multi-source abnormal information in Beidou system real-time orbit determination according to claim 1, characterized in that: In step S5, the detection factors are clustered by the DBSCAN clustering analysis method. The specific method for obtaining the abnormal detection factors includes the following sub-steps: S4-1, set the neighborhood radius R and the minimum number of points in the neighborhood minpoints; S4-2, traverse all the detection factors of the epoch to be detected, and for the detection factors whose number of points within the radius is greater than or equal to minpoints, include them as core points in the core point list; form a temporary cluster corresponding to the core point with the detection factors within the radius of the core point; S4-3, for each temporary cluster, check whether the detection factor therein is a core point. If so, merge the temporary cluster corresponding to the detection factor as the core point with the current temporary cluster, until each detection factor in the temporary cluster is not in the core point list, or the detection factors within the neighborhood radius R of each detection factor in the temporary cluster are in the same temporary cluster, and obtain a cluster; S4-4. Detection factors that are not classified into clusters are regarded as abnormal detection factors.

9. The method for real-time detection of multi-source abnormal information for Beidou system real-time orbit determination according to claim 1, characterized in that: The specific method of determining the source of the abnormality by calculating the proportion of the abnormality detection factor in step S6 includes the following sub-steps: S5-1. Calculate the proportion of the anomaly detection factor, and its expression is: in is the proportion of anomaly detection factors; represents the total number of detection factors; represents the total number of anomaly detection factors; S5-2, when there is a filtered innovation value greater than or equal to the phase jump detection threshold in the epoch to be detected, and the proportion of the abnormal detection factor is 100%, it is determined that the corresponding satellite has a phase jump fault; When there is no filtering innovation value greater than or equal to the phase jump detection threshold in the epoch to be detected, and the proportion of the abnormal detection factor of each satellite in the Beidou system is less than 90%, it is determined to be an observation anomaly; When there is no filtering innovation value greater than or equal to the phase jump detection threshold in the epoch to be detected, and the proportion of the satellite's abnormal detection factor is greater than or equal to 90% and less than 100%, the satellite is judged to be in an abnormal state, and the abnormality is determined to be caused by the satellite atomic clock frequency jump; When there is no filtered innovation value greater than or equal to the phase jump detection threshold in the epoch to be detected, and the satellite's abnormal detection factor accounts for 100%, the satellite is judged to be in an abnormal state, and the abnormality is determined to be caused by satellite orbit maneuvers.

Images