CN117370918B - High-rail safety-oriented space non-cooperative target maneuvering passive detection method - Google Patents
High-rail safety-oriented space non-cooperative target maneuvering passive detection method Download PDFInfo
- Publication number
- CN117370918B CN117370918B CN202311671226.XA CN202311671226A CN117370918B CN 117370918 B CN117370918 B CN 117370918B CN 202311671226 A CN202311671226 A CN 202311671226A CN 117370918 B CN117370918 B CN 117370918B
- Authority
- CN
- China
- Prior art keywords
- detection
- sequence
- matrix
- sliding window
- time
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000001514 detection method Methods 0.000 title claims abstract description 107
- 238000005259 measurement Methods 0.000 claims abstract description 41
- 238000011156 evaluation Methods 0.000 claims abstract description 24
- 230000003287 optical effect Effects 0.000 claims abstract description 21
- 230000000737 periodic effect Effects 0.000 claims abstract description 21
- 238000010183 spectrum analysis Methods 0.000 claims abstract description 7
- 239000011159 matrix material Substances 0.000 claims description 53
- 238000000034 method Methods 0.000 claims description 36
- 230000006870 function Effects 0.000 claims description 22
- 238000004422 calculation algorithm Methods 0.000 claims description 20
- 239000013598 vector Substances 0.000 claims description 20
- 238000000354 decomposition reaction Methods 0.000 claims description 19
- 230000002159 abnormal effect Effects 0.000 claims description 14
- 238000004364 calculation method Methods 0.000 claims description 11
- 230000008859 change Effects 0.000 claims description 11
- 230000001788 irregular Effects 0.000 claims description 8
- 238000012545 processing Methods 0.000 claims description 8
- 238000005070 sampling Methods 0.000 claims description 8
- 238000004458 analytical method Methods 0.000 claims description 6
- 238000010586 diagram Methods 0.000 claims description 6
- 230000005856 abnormality Effects 0.000 claims description 4
- 230000008569 process Effects 0.000 claims description 4
- 230000001360 synchronised effect Effects 0.000 claims description 4
- 238000010276 construction Methods 0.000 claims description 3
- 230000004069 differentiation Effects 0.000 claims description 3
- 239000000203 mixture Substances 0.000 claims description 3
- 238000010606 normalization Methods 0.000 claims description 3
- 238000006467 substitution reaction Methods 0.000 claims description 3
- 230000002123 temporal effect Effects 0.000 claims description 3
- 230000000007 visual effect Effects 0.000 claims description 2
- 230000008447 perception Effects 0.000 claims 1
- 238000000605 extraction Methods 0.000 abstract description 2
- 238000001914 filtration Methods 0.000 description 4
- 238000004088 simulation Methods 0.000 description 4
- 239000000284 extract Substances 0.000 description 2
- 238000004519 manufacturing process Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012216 screening Methods 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 230000009471 action Effects 0.000 description 1
- 230000001174 ascending effect Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 230000002950 deficient Effects 0.000 description 1
- 238000005265 energy consumption Methods 0.000 description 1
- 239000000446 fuel Substances 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000007689 inspection Methods 0.000 description 1
- 230000035772 mutation Effects 0.000 description 1
- 238000000513 principal component analysis Methods 0.000 description 1
- 230000000630 rising effect Effects 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
- G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
- G01C21/24—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 specially adapted for cosmonautical navigation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/10—Pre-processing; Data cleansing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
- G06F18/213—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods
- G06F18/2135—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods based on approximation criteria, e.g. principal component analysis
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/24—Classification techniques
- G06F18/243—Classification techniques relating to the number of classes
- G06F18/2433—Single-class perspective, e.g. one-against-all classification; Novelty detection; Outlier detection
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Computer Vision & Pattern Recognition (AREA)
- General Engineering & Computer Science (AREA)
- Mathematical Physics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Artificial Intelligence (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Radar, Positioning & Navigation (AREA)
- Evolutionary Biology (AREA)
- Evolutionary Computation (AREA)
- Remote Sensing (AREA)
- Mathematical Analysis (AREA)
- Computational Mathematics (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Astronomy & Astrophysics (AREA)
- Computing Systems (AREA)
- Automation & Control Theory (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- Navigation (AREA)
Abstract
The invention discloses a high-orbit safety-oriented space non-cooperative target maneuvering passive detection method, which is characterized in that a two-dimensional relative sight line measurement angle sequence is obtained based on an optical camera to serve as input, the periodic component of the detection sequence is fully utilized by means of the characteristic extraction capability of singular spectrum analysis, a differential evaluation index is constructed, a maneuvering threshold value is adaptively selected, the detection results of the periodic component, the half-period component and the two-dimensional angle measurement sequence are fused by an AND gate and an OR gate respectively, and maneuvering time points are comprehensively judged, so that maneuvering detection of any target on a GEO type orbit is realized. According to the invention, other priori information and assumption conditions are not needed, and the maneuvering detection judgment of the space-based non-cooperative target can be performed only by using the two-dimensional angle measurement sequence output by the space-based angle measurement system.
Description
Technical Field
The invention relates to the technical field of space-based passive detection, in particular to a space non-cooperative target maneuvering passive detection method for high-rail safety early warning.
Background
In the age background of space commercialization, which is increasingly intensified and numerous launching tasks, the cost for manufacturing the small satellites is reduced so that mass production can be carried out. Because of the limited radio frequency and orbital resources of satellites, most of the available and useful satellite frequency and orbital resources are already occupied, with the space resources of communication and navigation satellites being particularly scarce. Current world orbits and spectrum resource competition have tended to become white-hot, with the consequent increased spatial security problems. Space debris, invalid satellites, satellite weapons and other space non-cooperative targets have more and more threats to the on-orbit active spacecraft, and the exponential rise of maneuvering times of the space debris, invalid satellites, satellite weapons and other space non-cooperative targets has great threats to the safety of the spacecraft on-orbit operation. In this context, it is important for non-cooperative targets to have the ability to detect, sense, and take action on their orbital impulse maneuvers in time.
Currently, for the detection task of a space-based non-cooperative target, the relative orbit determination system based on passive angle measurement of an optical camera is widely applied. Compared with other commonly used spaceborne measurement sensors, the optical camera has the characteristics of high concealment, low energy consumption, simplicity, convenience, reliability and the like, and is particularly suitable for completing the passive detection task of a space non-cooperative target.
In an angle-only observation system on a space base, a conventional maneuvering detection method is generally applied to a tracking filtering stage, namely, an innovation sequence or distance information in the filtering method is used as detection quantity to screen abnormal points; however, in the initial orbit determination stage, the input observed quantity only has two angle values, no other priori information exists, and the characteristics of the satellite under different orbits and different initial relative orbits affect the angle observed quantity, namely, the change rule of the angle observed quantity is different, so that the maximum/small value, average value and other characteristic quantities cannot be used, and the maneuvering detection task in the stage is difficult. However, in a certain period of time, the satellites flying freely meet a certain orbit dynamics rule, so that the change trend of the angle observance quantity also meets a certain rule. Under the condition of no other priori information, the maneuvering detection method at the stage is constructed by focusing on whether the target maneuvers or not by analyzing the change trend of the angle observance.
With a general time sequence (a set of observations y 1 , …, y n And arranged in time sequence), the basic model of the angle observation quantity is an addition model, and has three changing characteristics in the following formula:
wherein the trending componentThe change trend of the angle measurement sequence under the influence of various factors in a long period is reflected; irregular component +.>Reflecting error sums such as observation errors of the angle measurement sequence; whereas anomalies caused by maneuver of the target satellite usually correspond to periodically varying features +.>Is a sudden change of the number of the (a). The mobility detection problem at this time is converted into a problem of extracting periodic feature mutations of the time series.
The current time sequence anomaly detection algorithm applied to the space-based angle measurement system comprises the following steps: the problems of target maneuver detection can be solved by principal component analysis, wavelet change, empirical mode decomposition and the like, but the anomaly detection methods have higher requirements on priori information, or require a large amount of historical information as input to finish accurate maneuver moment screening, or finish an initial orbit determination task, and utilize abundant orbit parameters as test statistics in a tracking filtering stage to realize maneuver detection. In the initial orbit determination stage in the angle observation of the space base, a group of observation quantity only has two angle values, no prior conditions such as a parameter model, a assumed stationarity condition and the like are assumed, and the application difficulty of the single anomaly detection algorithm is high.
Disclosure of Invention
Aiming at the problems in the existing anomaly detection algorithm, the invention provides a spatial non-cooperative target maneuvering passive detection method for high-rail safety early warning, which only uses angle data under the condition of insufficient prior inspection information and carries out continuous and accurate maneuvering detection and discrimination on a target in real time under the condition of light-weight satellite-borne calculation load.
A passive detection method for a spatial non-cooperative target maneuver facing high-rail safety comprises the following steps:
step 1, installing an optical camera at the centroid of a perceived satellite, and establishing an optical camera relative sight measurement model under a local vertical and local horizontal coordinate system to obtain a two-dimensional relative sight measurement angle sequenceAs a detection sequence; step 2, defining a singular spectrum analysis algorithm, and defining a proper safety interval and a sliding window length valueFor a two-dimensional detection sequence +.>Performing two-channel synchronous singular value decomposition in the current sliding window to obtain time sequence modeling of different trend components, namely real-time empirical orthogonal functions; step 3, defining a weight correlation coefficient calculation algorithm, and removing trending components, irregular components and the like and calculating periodic components and half-period components after reconstruction by calculating weight correlation coefficients among different time experience orthogonal functions to adaptively classify; step 4, defining a differential evaluation index, a basic window length and a hysteresis window length, respectively calculating the differential evaluation index of a period item and a half period item, setting a detection threshold value, and judging whether the change of the differential evaluation index is caused by a target maneuver or not based on the detection threshold value; and 5, deploying the angle measurement maneuvering detection algorithm on a sensing satellite, and inputting the relative measurement angle measured by the optical camera into the algorithm to realize on-line discrimination of maneuvering detection of the target satellite.
Preferably, the optical camera relative line of sight measurement model established in step 1 specifically includes:wherein->Representing the tracking spacecraft A and the target spacecraft B in LVLH system, respectively>A position vector of the moment.
Azimuth anglePitch angle->As inputs to the subsequent steps, they are synchronized independently of each other, and the integrated processing is performed only in step 5.
Preferably, the singular spectrum analysis algorithm in the step 2 does not need prior conditions such as a parameter model, a stability condition and the like, is suitable for an initial orbit determination stage with insufficient prior information, and extracts different components in a time sequence by performing operations such as decomposition, reconstruction and the like on a track matrix of the time sequence to be researched, so that tasks such as feature analysis, denoising and the like are completed.
Preferably, the step 2 specifically comprises: step 2.1, after the length of the safety interval is set, carrying out sliding window processing on a certain dimension angle sequence obtained in the step 1; step 2.2, embedding the detection sequences in each sliding window in sequence, and performing hysteresis arrangement to obtain a track matrix; and 2.3, carrying out singular value decomposition on the track matrix, and arranging the track matrix from large to small according to the singular values to obtain time experience orthogonal functions corresponding to different trend components.
Preferably, the sliding window length in step 2.1 is set to N. In step 2.2, the number of the track matrix lines is selected as L, and the detection sequences are arranged in a hysteresis way, so thatThen track matrix->Is->Is a matrix of (a):it is made of->Hysteresis vector of length L:composition is prepared.
Preferably, in step 2.3, in order to avoid the problems of different dimensions and larger calculation amount possibly encountered by singular value decomposition, an equivalent substitution method is selected, and a square matrix of the track matrix is calculated:and decomposing the characteristic value to obtain the characteristic value: />And corresponding feature vectors: />，/>Reflecting the evolution of the time series, called time empirical orthogonal functions.
Preferably, in step 3, correlation coefficients are calculated by using time empirical orthogonal functions of different trend components obtained by decomposition in pairs, and for a time sequence, the specific calculation method of the weight correlation coefficients is as follows: for time seriesWhereinIts weight correlation coefficient->. The cross-correlation coefficients of the L time empirical orthogonal functions are calculated and visualized-the weight correlation coefficient thermodynamic diagrams thereof are drawn.
Preferably, the weight correlation coefficient thermodynamic diagram is only used as a part of the visual analysis, and the specific selection process is as follows: calculating a hysteresis vectorAt->Projection on, i.e. the temporal principal components: />Reconstruction is performed by a time empirical orthogonal function and a time principal component:. Setting cross-correlation coefficient threshold +.>The cross-correlation coefficient satisfies->The two sequences are regarded as a group, the sequences are orderly arranged from the big to the small according to the characteristic value, and the number of each sequence group is marked as +.>RemovingHigh frequency trending component of +.>Will be +.>And->The sequences within the group are weighted and summed to sequentially correspond to the periodic characteristic component and the half-period characteristic component of the current sliding window sequence.
Preferably, in step 4, the defined differential evaluation index calculation flow is as follows: setting the length of the basic window asHysteresis window length +.>The number of sampling windows is +.>And (3) simultaneously carrying out the operations of the steps 2.2-2.3 on the detection sequences in the two windows to respectively obtain:
basic track matrixAnd its characteristic valueAnd corresponding feature vector->；
Hysteresis trajectory matrixAnd its characteristic valueCorresponding feature vector->；
Before selectingConstructing a differential evaluation index by the feature vectors: firstly, the hysteresis track matrix is compared with the basic track matrix, namely the +.>Is subjected to normalization processing:further, different weights are calculated according to singular values of the hysteresis track matrix at the same time, and a final differential evaluation index is solved in a weighted mode>。
Preferably, the specific procedure for adaptively selecting the detection threshold is as follows: calculating covariance among row vectors of track matrix in current sliding window, and constructing auto-covariance matrix，/>The main diagonal element is the variance of the current window sequence. According to->Reconstructing a main component +.>And mean value normalized by irregular component ++>. Threshold value of the method for constructing maneuver detection>：/>Wherein->Is composed of->And calculating according to a differential construction index calculation method. When->When the detection sequence in the current sliding window is judged to be abnormal, namely the current moment is judged to be an abnormal point, and the logic calibration is 1; on the contrary, when->When the logic is calibrated to 0.
Preferably, the periodic characteristic component and the periodic characteristic component obtained in the step 3 are independently synchronized as the input of the step 4, the judgment result of the periodic characteristic component and the periodic characteristic component is subjected to AND gate logic judgment, and the final result is output: 0 indicates that the current time is an abnormal point-free bit, and 1 indicates that the current time is an abnormal point.
Preferably, in step 5, two goniometric sequences (yaw anglePitch angle->) Respectively taking the two results as the input of the step 2, respectively and independently performing the operations of the step 2 to the step 4, and performing OR gate logic judgment on the judgment results of the two results, wherein each timeAnd outputting a final judging result corresponding to the logic of one-bit calibration: logic 0 indicates no abnormality in the current sliding window, and the central moment point is the moment point at which no maneuver occurs; logic 1 indicates that the data anomaly in the current sliding window is caused by maneuver, and the central moment point is the maneuver moment point.
The beneficial effects are that:
(1) The invention solves the problem that prior information in the space-based passive angle measurement data is deficient and the maneuvering detection in the relative orbit determination process can not be realized, utilizes the characteristics of low prior information requirement, strong data decomposition capability and the like in the singular spectrum analysis to perform characteristic analysis on the data, and quantizes the periodic component of the detection sequence, thereby realizing real-time maneuvering detection on the target by utilizing the differential evaluation index;
(2) The invention selects the maneuvering detection threshold value in a self-adaptive way instead of manual selection, thereby effectively improving the success rate of maneuvering detection and having higher applicability in engineering practice;
(3) In the invention, in the step 4, the periodic characteristic components and the half-period characteristic components are fully utilized, an AND gate logic discriminator is added, and the comprehensive judgment is carried out to reduce the false alarm rate; in the step 5, the detection results of the two angle measurement sequences are comprehensively processed, an OR gate logic discriminator is added to reduce the false alarm rate, and meanwhile, the universality of the application of the maneuvering detection algorithm on various tracks is improved;
(4) The method fully improves the algorithm to improve the accuracy of the maneuvering detection of the target without introducing extra equipment;
(5) The invention is applicable to GEO type rails, has wider freedom degree for rail guidance and does not need to increase fuel cost;
(6) The invention has the advantages of no need of changing the traditional structure of the satellite sensor, no limitation to the relative distance of the satellite, and wider freedom of maneuvering detection of relative targets.
Drawings
FIG. 1 is a schematic diagram of a space-based optical camera measurement geometry according to one embodiment of the present invention;
FIG. 2 is a logic flow diagram of a maneuver detection algorithm according to an embodiment of the present invention;
FIGS. 3 (a) and 3 (b) are graphs of window changes of the differential evaluation indexes 25s and 10s under different track matrix window lengths according to one embodiment of the present invention;
FIG. 4 is a graph of detection success rate and detection error at different relative track heights and sensor accuracy for one embodiment of the present invention;
FIG. 5 is a graph of detection success rate and detection error at different maneuver levels and sensor accuracy for one embodiment of the present invention.
FIG. 6 is a graph showing comparison of detection success rate and detection error of empirical mode decomposition methods at different maneuver levels and methods of the present invention, according to an embodiment of the present invention.
Detailed Description
The following description of the embodiments of the present invention will be made clearly and completely with reference to the accompanying drawings, in which it is apparent that the embodiments described are only some embodiments of the present invention, but not all embodiments. All other embodiments, which can be made by those skilled in the art based on the embodiments of the invention without making any inventive effort, are intended to be within the scope of the invention.
The invention discloses a space non-cooperative target maneuvering passive detection method for high-orbit safety precaution, aiming at the problems that the prior space-based passive maneuvering detection algorithm has higher requirement on priori information, a large amount of historical information is needed as input to finish accurate maneuvering time screening, the maneuvering detection is usually realized by using rich orbit parameters as test statistics in a tracking filtering stage, but in an initial orbit determination stage, a group of observed quantities in space-based angle observation only have two angle values, no priori conditions such as a parameter model, a assumed stability condition and the like are assumed, and the application difficulty of the single anomaly detection algorithm is higher; the periodic components of the detection sequences are fully utilized by means of the characteristic extraction capability of singular spectrum analysis, the differential evaluation indexes are built, the maneuvering threshold value is selected in a self-adaptive mode, the detection results of the periodic and half-period components and the two-dimensional angle measurement sequences are fused by means of an AND gate and an OR gate respectively, maneuvering time points are comprehensively judged, and therefore maneuvering detection of any target on a GEO type track is achieved. According to the invention, other priori information and assumption conditions are not needed, and the maneuvering detection judgment of the space-based non-cooperative target can be performed only by using the two-dimensional angle measurement sequence output by the space-based angle measurement system.
In an initial orbit determination stage, a spatial non-cooperative target maneuvering passive detection method facing high-orbit safety early warning is provided, and the method comprises the following steps:
step 1, installing an optical camera at the centroid of a perceived satellite, and establishing an optical camera relative sight measurement model under a local vertical and local horizontal coordinate system to obtain a two-dimensional relative sight measurement angle sequenceAs a sequence of detection amounts.
Step 1.1, an optical camera relative vision measurement model shown in fig. 1 is established, and the expression is specifically as follows:wherein->Representing the tracking spacecraft A and the target spacecraft B in LVLH system, respectively>A position vector of the moment.
Step 1.2, as shown in FIG. 2, the azimuth angles are synchronously input independently of each otherPitch angle->In step 5, the integrated process is performed.
Step 2, constructing a singular spectrum analysis algorithm: the method does not need prior conditions such as a parameter model, a stability condition and the like, is suitable for an initial orbit determination stage with insufficient prior information, and extracts different components in a time sequence by performing operations such as decomposition, reconstruction and the like on a track matrix of the time sequence to be researched, thereby completing any of feature analysis, denoising and the likeAnd (5) carrying out business. Defining proper safety interval and sliding window length value, and for two-dimensional detection quantity sequenceAnd carrying out the singular value decomposition of the double-channel synchronization in the current sliding window to obtain time sequence variants of different trend components, namely, real-time empirical orthogonal functions.
And 2.1, after the set safety interval length is 100s, performing sliding window processing on the certain dimensional angle sequence obtained in the step 1, wherein the sliding window length is set to be N=100deg.S.
Step 2.2, embedding the detection sequences in each sliding window in sequence, and obtaining a track matrix by hysteresis arrangement: selecting the length of a track matrix window as L, and performing hysteresis arrangement on the detection sequence to enableThen track matrix->Is->Is a matrix of (a): />It is made of->Hysteresis vector of length L:composition is prepared.
Step 2.3, performing singular value decomposition on the track matrix, and arranging the track matrix from large to small according to the singular values to obtain time experience orthogonal functions corresponding to different trend components: in order to avoid the problems of different dimensions and larger calculated amount possibly encountered by singular value decomposition, an equivalent substitution method is selected, and a square matrix of a track matrix is calculated:decomposing the characteristic value to obtain the characteristic value：/>And corresponding feature vectors: />，/>Reflecting the evolution of the time series, called time empirical orthogonal functions.
Step 3, defining a weight correlation coefficient calculation algorithm, and removing trending components, irregular components and the like and calculating periodic components and half-period components after reconstruction by calculating weight correlation coefficients among different time experience orthogonal functions to adaptively classify, wherein the method specifically comprises the following steps of:
step 3.1, calculating correlation coefficients of time experience orthogonal functions of different trend components obtained through decomposition in pairs, wherein for a time sequence, the specific calculation method of the weight correlation coefficients comprises the following steps: for time seriesIts weight correlation coefficient->. The cross-correlation coefficients of the L time-empirical orthogonal functions are calculated.
Step 3.2, selecting periodic components: calculating a hysteresis vectorAt->Projection on, i.e. the temporal principal components:reconstruction is performed by a time empirical orthogonal function and a time principal component: />。
Step 4, defining a differential evaluation index, a basic window length and a hysteresis window length, respectively calculating the differential evaluation index of a period item and a half period item, setting a detection threshold value, and judging whether the change of the differential evaluation index is caused by target maneuver or not based on the detection threshold value, wherein the specific steps are as follows:
step 4.1, a defined differential evaluation index calculation flow is as follows: setting the length of the basic window asHysteresis window length +.>The number of sampling windows is +.>And (3) simultaneously carrying out the operations of the steps 2.2-2.3 on the detection sequences in the two windows to respectively obtain:
basic track matrixAnd its characteristic valueAnd corresponding feature vector->。
Hysteresis trajectory matrixAnd its characteristic valueCorresponding feature vector->。
Step 4.2, before selectionConstructing a differential evaluation index by the feature vectors: firstly, the hysteresis track matrix is compared with the basic track matrix, namely the +.>Is subjected to normalization processing:further, different weights are calculated according to singular values of the hysteresis track matrix at the same time, and a final differential evaluation index is solved in a weighted mode>。
Step 4.3, defining a self-adaptive detection threshold selecting method, wherein the specific flow is as follows: calculating covariance among row vectors of track matrix in current sliding window, and constructing auto-covariance matrix，/>The main diagonal element is the variance of the current window sequence. According to->Reconstructing principal components that substantially reflect the overall characteristics of the original sequence from partial eigenvalues and eigenvectors of (a)And mean value normalized by irregular component ++>. Threshold value of the method for constructing maneuver detection>：/>Wherein->Is composed of->Press differenceAnd calculating by using a dissimilarity construction index calculation method. When->When the detection sequence in the current sliding window is judged to be abnormal, namely the current moment is judged to be an abnormal point, and the logic calibration is 1; conversely, whenWhen the logic is calibrated to 0.
When comprehensively judging whether the single angle measurement sequence is abnormal, independently synchronizing the periodic characteristic component and the half-period characteristic component obtained in the step 3 as the input of the step 4, performing AND gate logic judgment on the judgment results of the periodic characteristic component and the half-period characteristic component, and outputting a final result: 0 indicates that the current time is an abnormal point-free bit, and 1 indicates that the current time is an abnormal point.
Step 5, two goniometric sequences (yaw anglePitch angle->) And (3) respectively taking the judgment results of the step (2) and the step (4) as the input of the step (2), respectively and independently carrying out OR gate logic judgment on the judgment results of the step (2) and the step (4), and outputting a final judgment result corresponding to one-bit calibration logic at each moment: logic 0 indicates no abnormality in the current sliding window, and the central moment point is the moment point at which no maneuver occurs; logic 1 indicates that the data anomaly in the current sliding window is caused by maneuver, and the central moment point is the maneuver moment point.
The angle measurement maneuvering detection algorithm is deployed on a sensing satellite, each group of relative measurement angles with certain continuous time measured by an optical camera is input into the algorithm, a threshold value is selected in a self-adaptive mode, and multidimensional data are fused and judged, so that online judgment of maneuvering detection of a target satellite is realized.
The applicability of the present invention will be described with reference to the following examples.
The following calculation conditions and technical parameters are set:
1) The orbit semi-long axis of the perceived satellite A is 42166.3 km, the eccentricity is 0.0005, the orbit inclination angle is 0 degrees, the perigee amplitude angle is 0 degrees, the right-hand warp of the ascending intersection point is 0 degrees, and the true-nearby point angle is 0 degrees;
2) The orbit semi-long axis of the target satellite is 42066.3 km, the eccentricity is 0.0005, the orbit inclination angle is 0/0.5 degrees, the near-site amplitude angle is 0 degrees, the right-hand warp of the rising intersection point is 0 degrees, and the true near-point angle is 0 degrees;
3) The mean square error of the angle measurement noise of the camera is 0.0005rad, and the mean square error of the absolute position noise of the sensing satellite is 100m;
4) The momentum of the target satellite is 0.5/0.75/1.0/1.25/1.5/1.75/2m/s respectively;
5) The relative distance between the target satellite and the perceived satellite is 500/625/750/875/1000km respectively;
6) The continuous time of the measurement angle sampling is 900s;
7) The track matrix window length L is 10s/25s respectively;
8) And when the success rate and the detection error are counted, the number of target simulation times of each track type is 500.
Based on the angle-measuring-only relative orbit determination method and the set calculation conditions and technical parameters, simulation verification is carried out.
As shown in fig. 3 (a) and 3 (b), the relative measurement angle sampling continuous time is 900s, the target satellite orbit inclination angle target is 0 degree, the momentum of the satellite machine is fixed to be 1m/s, the relative orbit height difference of the target satellites is 1000km, and the average square error of the angle measurement noise of the optical camera is 0.00003 rad time difference differentiation evaluation indexIs a graph of the variation of (a). As shown in the figure, when selecting the number of rows of different track matrixes, the differential evaluation index is +.>The influence of the variation of (c) is large. FIG. 3 (a) is the result of the trace matrix window length being 25s; FIG. 3 (b) shows the result of the track matrix window length being 10s, it is evident that the maneuver-induced abrupt change is more pronounced when the track matrix window length is smaller, but too small a window length results in too much computation, soThe selection of the threshold value is more convenient, so the length L=10s of the track matrix window is selected in the simulation experiment.
As shown in fig. 4, the relative measurement angle sampling continuous time is 900s, the target satellite orbit dip angle target is 0 °, the satellite machine momentum size is fixed to be 1m/s, the relative orbit height difference of the target satellites is 500/625/750/875/1000km, and the detection success rate and detection error curve of the angle measurement noise mean square error of the optical camera are 0.003/0.0005/0.0003/0.00003 rad respectively; as shown in FIG. 5, the relative measurement angle sampling continuous time is 900s, the target satellite orbit inclination angle target is 0.5 degrees, the target satellite relative orbit height difference is fixed to 1000km, the maneuvering amount of the target satellite is 0.5/0.75/1.0/1.25/1.5/1.75/2m/s, and the detection success rate and detection error curve of the angle measurement noise mean square error of the optical camera are 0.003/0.0005/0.0003/0.00003 rad.
As can be seen from the graphs in FIGS. 4-5, the detection accuracy of the non-cooperative target maneuver detection is higher by the method. Under the method, the detection success rate is more than 90% and the detection error is less than 80s under the condition of the dynamic quantity of the magnitude of 1m/s and the angular noise mean square error of the 0.0003 rad optical camera; after the track inclination angle is changed, 100% of detection is successful under the conditions of the maneuvering quantity of the order of 2m/s and the average square error of the angle measurement noise of the 0.0003 rad optical camera, and the detection error is within 30 s.
As shown in FIG. 6, the ratio of the detection success rate and the detection error curve of the method of the present invention are compared with the empirical mode decomposition method (Empirical Mode Decomposition, EMD) with the relative measurement angle sampling continuous time of 900s, the target satellite orbit inclination angle fixed at 0 degrees, the angular noise mean square error fixed at 0.0003 rad, the target satellite relative orbit height difference fixed at 625km, and the satellite momentum of 0.5/0.75/1.0/1.25/1.5/1.75/2 m/s. The EMD method detection data and the target simulation parameters are identical to the application method of the invention. As can be seen from fig. 6, in the same track height difference (set to 625 km), the detection success rate of the method of the present invention is improved by more than about 15% compared with the EMD method, and the detection error is reduced by about 50%, which indicates that the mobility detection capability of the method is significantly better than the EMD method in different mobility levels.
Therefore, by adopting the method, the maneuvering detection of the initial orbit determination stage of the non-cooperative target under the condition of no priori information can be realized by only depending on the passive angle measurement of the satellite-borne optical camera. In particular, compared with the traditional maneuvering detection method (such as an EMD method), the method can realize accurate detection of non-cooperative targets with various relative track types and different maneuvering levels under the condition of lacking prior information and ranging information, and is a breakthrough progress brought by the method.
Finally, it should be noted that: the foregoing description is only a preferred embodiment of the present invention, and the present invention is not limited thereto, but it is to be understood that modifications and equivalents of some of the technical features described in the foregoing embodiments may be made by those skilled in the art, although the present invention has been described in detail with reference to the foregoing embodiments. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention should be included in the protection scope of the present invention.
Claims (4)
1. The passive detection method for the spatial non-cooperative target maneuver facing the high-orbit safety is characterized by comprising the following steps of:
step 1, installing an optical camera at the centroid of a perceived satellite, and establishing an optical camera relative line-of-sight measurement model under a local vertical and local horizontal coordinate system:
wherein,respectively represent the t of the tracking spacecraft A and the target spacecraft B in the LVLH system i Position vector of time, n 1 、n 2 In order to measure the noise of the light,
obtaining a two-dimensional relative sight line measurement angle sequence { az (i), el (i) } as a detection quantity sequence to be detected respectively;
step 2, defining a singular spectrum analysis algorithm, carrying out two-channel synchronous singular value decomposition and reconstruction operation on a track matrix of a sequence to be researched in a current sliding window, extracting different components in a time sequence, carrying out feature analysis and denoising, and obtaining time sequence modeling types of different trend components, namely an empirical orthogonal function;
step 3, defining weight correlation coefficient calculation algorithm, self-adaptively classifying by calculating weight correlation coefficients among different time experience orthogonal functions, removing trending components and irregular components, calculating periodic components and half-period components after reconstruction, in particular,
calculating correlation coefficients of time experience orthogonal functions of different trend components obtained through decomposition in pairs: for time series x (i), y (i), where i=1, …, N is the sliding window length, the weight correlation coefficient is expressed as:in (1) the->Corresponding data average values in the selected window; calculating and visualizing the cross-correlation coefficients of the L time experience orthogonal functions, and drawing a weight correlation coefficient thermodynamic diagram of the L time experience orthogonal functions, wherein the weight correlation coefficient thermodynamic diagram is used as a part of visual analysis, and the specific selection process comprises the following steps: calculating a hysteresis vector X at the moment i i In U i Projection on, i.e. the temporal principal components:
reconstruction is performed by a time empirical orthogonal function and a time principal component: />Wherein a is m As a time principal component, U is a feature vector, subscript i represents track matrix row number, j represents track matrix column number, and m represents trendOrdinal number of the composition; setting a cross-correlation coefficient threshold kappa, wherein when the cross-correlation coefficient meets r & gtkappa, two sequences are regarded as a group, the two sequences are sequentially arranged from large to small according to characteristic values, the number of each sequence group is recorded as g, a high-frequency trend component with g=1 and an irregular component with g & gt4 are removed, and sequences in the groups with g=2 and g=3 are respectively weighted and summed to sequentially correspond to a periodic characteristic component and a half-period characteristic component of a current sliding window sequence;
step 4, constructing a differential evaluation index, respectively calculating the differential evaluation index of the period item and the half period item, setting a detection threshold value, judging whether the change of the differential evaluation index is caused by the target maneuver or not, specifically,
step 4.1, constructing a differentiation evaluation index: setting the length of the basic window asHysteresis window length +.>The number of sampling windows is delta, the detection sequences in the two windows are simultaneously subjected to the singular value decomposition operation in the step 2, and a basic track matrix F, a hysteresis track matrix H and a corresponding feature vector U are respectively obtained f 、U h The method comprises the steps of carrying out a first treatment on the surface of the The first eta feature vectors are selected to construct a differentiation evaluation index: solving U ht T=1,.. single change characteristic of eta and normalization processing: />Meanwhile, calculating different weights according to singular values of the hysteresis track matrix, and weighting and solving a final differential evaluation index Ne;
step 4.2, adaptively selecting a detection threshold value: calculating covariance among row vectors of track matrix in current sliding window, and constructing auto-covariance matrix T x The method comprises the steps of carrying out a first treatment on the surface of the According to T x Is used for calculating the principal component T of the integral features of the original sequence by partial eigenvalues and eigenvectors main And the normalized mean value of the irregular componentConstructing a threshold gamma of the maneuver detection method: />Wherein Ne γ Is composed of T main Calculated according to a differential construction index calculation method; when Ne is more than or equal to gamma, judging that the detection sequence in the current sliding window is abnormal, namely judging that the current moment is an abnormal point, and logically calibrating to be 1; conversely, when Ne < gamma, the logic is calibrated to be 0;
and 4.3, taking the periodic characteristic component and the periodic characteristic component obtained in the step 3 as the input of the step 4 respectively, performing AND gate logic judgment on two judgment results correspondingly output in the step 4, and outputting a final result: 0 represents that the current moment is an abnormal point-free bit, and 1 represents that the current moment is an abnormal point;
step 5, deploying an angle measurement maneuver detection algorithm only on a perception satellite, and measuring two angle measurement sequences by an optical camera: the yaw angle az and the pitch angle el are respectively used as the input of the step 2 to correspondingly obtain two judgment results of the step 4 for OR gate logic judgment, each moment corresponds to one-bit calibration logic, and a final judgment result is output: logic 0 indicates no abnormality in the current sliding window, and the central moment point is the moment point at which no maneuver occurs; logic 1 indicates that the abnormality of the data in the current sliding window is caused by maneuver, and the central moment point is the maneuver moment point, so that the online judgment of the maneuver detection of the target satellite is realized.
2. The high-rail-security-oriented spatial non-cooperative target maneuvering passive detection method according to claim 1, wherein the step 2 is specifically: step 2.1, setting a safety interval length and a sliding window length value, performing sliding window processing on the certain dimensional angle sequence obtained in the step 1, and setting two channels to be performed simultaneously so as to realize two-dimensional sequence processing; step 2.2, embedding the detection sequences in each sliding window in sequence, and performing hysteresis arrangement to obtain a track matrix; and 2.3, carrying out singular value decomposition on the track matrix, and arranging the track matrix from large to small according to the singular values to obtain time experience orthogonal functions corresponding to different trend components.
3. The method for mechanically and passively detecting a space non-cooperative target for high-track safety according to claim 2, wherein in step 2.2, the number of rows L of the track matrix is selected, the detection sequence is arranged in a hysteresis manner, so that k=n-l+1, N is a sliding window length, and the track matrix X is represented as a matrix of lxk:
i.e. the trajectory matrix X comprises K hysteresis vectors of length L:
X i ＝(x i ,...,x i+L-1 ) T ,1≤i≤K。
4. the high-rail-security-oriented spatial non-cooperative target maneuver-passive detection method as defined in claim 3, wherein in step 2.3, the square matrix of the trajectory matrix is calculated based on the equivalent substitution method: s=xx T And decomposing the characteristic value to obtain the characteristic value: lambda (lambda) 1 ＞λ 2 ＞…＞λ L 0 and corresponding feature vector: u= [ U ] 1 ,U 2 ,…,U L ]U reflects the evolution of the time series, called time-empirical orthogonal functions.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311671226.XA CN117370918B (en) | 2023-12-07 | 2023-12-07 | High-rail safety-oriented space non-cooperative target maneuvering passive detection method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311671226.XA CN117370918B (en) | 2023-12-07 | 2023-12-07 | High-rail safety-oriented space non-cooperative target maneuvering passive detection method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN117370918A CN117370918A (en) | 2024-01-09 |
| CN117370918B true CN117370918B (en) | 2024-04-12 |
Family
ID=89391393
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202311671226.XA Active CN117370918B (en) | 2023-12-07 | 2023-12-07 | High-rail safety-oriented space non-cooperative target maneuvering passive detection method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117370918B (en) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN118243112A (en) * | 2024-05-24 | 2024-06-25 | 南京航空航天大学 | Space-based passive cooperative multi-target positioning method in complex space environment |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108519110A (en) * | 2018-04-26 | 2018-09-11 | 北京空间飞行器总体设计部 | Space non-cooperative target Autonomous Relative Navigation validating in orbit system based on image information |
| CN109307862A (en) * | 2018-07-05 | 2019-02-05 | 西安电子科技大学 | A kind of target radiation source individual discrimination method |
| CN109375279A (en) * | 2018-10-16 | 2019-02-22 | 淮海工学院 | A kind of static weight observation data gravity tide correction extracting method |
| CN116698049A (en) * | 2023-08-08 | 2023-09-05 | 南京航空航天大学 | Maneuvering detection method for passive detection of initial orbit determination section of space non-cooperative target |
-
2023
- 2023-12-07 CN CN202311671226.XA patent/CN117370918B/en active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108519110A (en) * | 2018-04-26 | 2018-09-11 | 北京空间飞行器总体设计部 | Space non-cooperative target Autonomous Relative Navigation validating in orbit system based on image information |
| CN109307862A (en) * | 2018-07-05 | 2019-02-05 | 西安电子科技大学 | A kind of target radiation source individual discrimination method |
| CN109375279A (en) * | 2018-10-16 | 2019-02-22 | 淮海工学院 | A kind of static weight observation data gravity tide correction extracting method |
| CN116698049A (en) * | 2023-08-08 | 2023-09-05 | 南京航空航天大学 | Maneuvering detection method for passive detection of initial orbit determination section of space non-cooperative target |
Also Published As
| Publication number | Publication date |
|---|---|
| CN117370918A (en) | 2024-01-09 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106372646B (en) | Multi-target tracking method based on SRCK-GMCPHD filtering | |
| CN105371870B (en) | A kind of in-orbit accuracy measurement method of star sensor based on star chart data | |
| CN104778358B (en) | The partly overlapping extension method for tracking target in monitored area be present in multisensor | |
| CN117370918B (en) | High-rail safety-oriented space non-cooperative target maneuvering passive detection method | |
| CN116450711B (en) | GNSS coordinate time sequence data stream matching method | |
| CN111798491A (en) | Maneuvering target tracking method based on Elman neural network | |
| CN103247057A (en) | Road target multi-hypothesis tracking algorithm under target-echo-road network data association | |
| CN105021199B (en) | Multi-model self-adapting method for estimating state and system based on LS | |
| Kong et al. | Application of stabilized numerical integration method in acceleration sensor data processing | |
| CN113411744A (en) | High-precision indoor positioning and tracking method | |
| CN115523927A (en) | GEO spacecraft maneuvering detection method based on optical sensor observation | |
| CN111830501A (en) | HRRP (high resolution representation protocol) historical characteristic assisted signal fuzzy data association method and system | |
| CN106501815B (en) | Space target track maneuvering fusion detection method only based on space-based angle measurement tracking | |
| CN113030940B (en) | Multi-star convex type extended target tracking method under turning maneuver | |
| CN115343744A (en) | Optical single-double-star combined on-satellite positioning method and system for aerial moving target | |
| CN115372957A (en) | Trajectory tracking method for hypersonic aircraft | |
| CN117271979A (en) | Deep learning-based equatorial Indian ocean surface ocean current velocity prediction method | |
| CN115963503A (en) | Method for quickly measuring relative pose of space non-cooperative slow-rotation target by laser radar | |
| CN116026325A (en) | Navigation method and related device based on neural process and Kalman filtering | |
| CN116108314A (en) | Method for determining space target track based on multi-sensor combination measurement | |
| Zhu et al. | An Accurate Line-of-Sight Rate Estimation Method Based on LSTM Recurrent Neural Network for Strapdown Imaging Seeker | |
| Muratov et al. | Use of AI for Satellite Model Determination from Low Resolution 2D Images | |
| Wang et al. | Statistical confidence domain data driven based fast in-flight alignment method | |
| CN118246252B (en) | Pulse maneuver detection method for non-cooperative target satellite | |
| CN110889227B (en) | Aircraft fuel measurement method based on multi-sensor information fusion |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |