Space target transaction discrimination method based on imaging auxiliary data

Abstract

The invention discloses a space target transaction discrimination method based on imaging auxiliary data, and the method comprises the steps: firstly obtaining the theoretical and measurement data of the azimuth angle and pitch angle of a space target rotary table, and eliminating outliers through employing an improved amplitude limiting filtering algorithm; filtering estimation is carried out on the data after the outliers are removed, and estimated real data are obtained to correct the measurement data; calculating and estimating the actual similarity between the real data and the theoretical data, and fitting a similarity function; for continuous new moments, the absolute value of the difference between the actual similarity and the estimated similarity obtained through the similarity function is calculated as a similarity difference value. According to the method, astronomical positioning is not needed, the target orbit transaction can be rapidly detected while target characteristic imaging is carried out, the research and judgment timeliness and the space target in-orbit operation safety are improved, the space target collision risk is effectively reduced, and the blank that the in-orbit satellite small-caliber telescope is used for detecting the target orbit transaction is filled.

Description

Technical Field

[0001] The present invention relates to the technical field of space situation awareness and space target recognition, and in particular to a method for distinguishing abnormal movement of space targets based on imaging auxiliary data. Background Art

[0002] With the rapid development of space technology, the scope of human activities has expanded from the traditional land, sea and air to space. Due to its special location advantage, space has become the core of communication, economy, military strategy and international relations on the earth. Now, it is the latest stage for human exploration, development and possible conquest. With the continuous increase of human activities in the space field and the continuous development and utilization of space resources, the number of space targets is growing exponentially.

[0003] In terms of low-orbit, with the construction and deployment of a series of giant constellations, space will be severely squeezed. Such a large number of densely distributed satellites will bring huge collision risks to the launch, normal operation and deorbit disposal of spacecraft. In the early stage of launch into orbit, Starlink satellites will be lifted for nearly one month under the action of small thrust, and there will be a phenomenon of thrust value jump before and after starting thrust and stopping thrust. At the same time, in order to avoid the risk of collision within the constellation, Starlink satellites are equipped with autonomous collision avoidance systems, which can receive space target cataloging data from the US Department of Defense and perform orbital maneuvers when necessary. Starlink satellites often make many small-scale orbital corrections and adjustments to avoid collisions, but Space Exploration Company does not fully release the situation of Starlink satellites performing orbital maneuvers, which will multiply the operating risks of other spacecraft. In terms of high-orbit, as satellite platforms develop from large-scale integration to small-scale decentralization, high-orbit satellites perform orbital maneuvers more frequently to carry out tasks such as satellite life extension and space experiments.

[0004] At present, spacecrafts of various countries frequently maneuver in orbit, especially when there is a collision warning event between domestic and foreign spacecrafts, the two sides often do not communicate, and both maneuver in a "back-to-back" manner, which poses a risk of collision due to maneuvers. The observation, monitoring and identification of these orbital anomaly space targets has always been the focus of research. Timely discovery and identification of orbital anomalies of such targets can effectively support space collision warning. The small-aperture telescopes carried by my country's in-orbit spacecraft can track targets for a long time and obtain a large amount of observation data. However, due to the equipment capabilities, the detection capabilities of such small-aperture telescopes are insufficient, resulting in a limited number of stars in the field of view. During the target tracking process, astronomical positioning measurements are often unable to be completed. At the same time, the angular measurement accuracy of such telescopes is poor, and it is difficult to accurately locate. Therefore, it is often impossible to use traditional anomaly detection methods for anomaly identification, and the data is difficult to fully utilize. Therefore, an effective method for identifying anomalies of space targets is urgently needed.

[0005] The principle of traditional motion detection for space targets is to compare the target's expected angular position in the camera coordinate system with the actual angular position, and consider motion to have occurred when it exceeds a certain detection threshold. Based on the target's predicted position and the coordinate system information of the monitoring satellite / camera, the Whether the target is visible relative to the camera at the time, and if it is visible, the position of the target in the camera coordinate system. The difference between the actual phase plane position and the predicted phase plane position of the target, in degrees, is recorded as , whether the target appears in the camera's field of view after the maneuver, assuming that the target is visible, depends on the size of the field of view, the difference between the predicted position and the actual position, and the direction of the predicted position. If the predicted position is at the edge of the field of view, when Smaller than the field of view When the position difference direction is ideal, the target will appear in the field of view after maneuvering; if the predicted position is in the center of the field of view, when Smaller than the field of view , the target will appear in the field of view after maneuvering.

[0006] Traditional abnormality detection uses feature-based The method is that the measured value and the theoretical predicted value are converted and calculated according to the following formula to obtain the characteristic quantity. If the value exceeds the set abnormality threshold, it is judged as an abnormality. Therefore, the main influencing factors of traditional abnormality detection and identification include the characteristic quantity observation error and the characteristic quantity theoretical calculation error. Specifically, the error is calculated using the astronomical positioning measured value and the astronomical positioning theoretical value that reflect the characteristic quantity.

[0007]

[0008] in, represents the measured value of astronomical positioning right ascension, It represents the theoretical value of astronomical positioning right ascension. represents the measured value of astronomical positioning declination, Indicates the theoretical value of astronomical positioning declination.

[0009] Although this method has the advantage of high measurement accuracy, its disadvantages are that the method is complicated and has poor timeliness. It needs to complete astronomical positioning first, requires high angle measurement accuracy, and is not suitable for imaging equipment. Summary of the invention

[0010] The purpose of the present invention is to overcome the defects and shortcomings of the prior art and to provide a method for distinguishing space target movement based on imaging auxiliary data.

[0011] To achieve the above object, the present invention provides the following technical solutions: A method for distinguishing space target movement based on imaging auxiliary data, characterized in that it specifically comprises the following steps: S1. At each moment, calculate the theoretical data of the turntable azimuth angle and the theoretical data of the turntable pitch angle of the space target according to the target orbit, and obtain the turntable azimuth angle measurement data and the turntable pitch angle measurement data of the space target during the tracking and imaging process; S2. The obtained turntable azimuth angle measurement data and turntable pitch angle measurement data are respectively subjected to outlier elimination using an improved limiting filter algorithm; S3. Perform filtering estimation on the turntable azimuth angle measurement data and the turntable pitch angle measurement data after outliers are eliminated, respectively, to obtain the turntable azimuth angle estimation true data and the turntable pitch angle estimation true data corresponding to each moment; S4. Using the obtained theoretical data of the turntable azimuth angle, the estimated real data of the turntable azimuth angle, the theoretical data of the turntable pitch angle and the estimated real data of the turntable pitch angle at each time, calculate the actual similarity between the estimated real data and the theoretical data at the corresponding time; S5. According to the calculated actual similarities of all moments, a similarity function is fitted to calculate the estimated similarity of the moment; S6. for each of the N consecutive new moments, obtain the actual similarity of the moment, use the similarity function to obtain the estimated similarity of the moment, and calculate the absolute value of the difference between the actual similarity and the estimated similarity at the moment as the similarity difference value, where N is a natural number greater than 0; S7. Determine whether the similarity difference values ​​of the N consecutive new moments meet the change condition, if so, determine that the space target has changed; if not, return to step S6.

[0012] Furthermore, the step S1 specifically includes: using the orbit of the space target and the azimuth and pitch measurement values ​​of the turntable to calculate the theoretical azimuth angle data and the theoretical pitch angle data of the turntable of the space target at each moment, and obtaining the turntable azimuth angle measurement data of the space target at each moment output in real time by the turntable azimuth encoder, and the turntable pitch angle measurement data of the space target at each moment output in real time by the turntable pitch encoder.

[0013] Furthermore, the step S2 specifically includes the following steps: S21. Based on the limiting filter algorithm, the following initial formula for outlier removal is proposed: (1) In formula (1), represents the data after outliers are removed. Indicates The sampling data corresponding to the moment, Indicates The sampling data corresponding to the moment, is a natural number greater than or equal to 1; Indicates the threshold value corresponding to the sampling data. When the sampling data is the turntable azimuth measurement data, the data after outlier removal is the turntable azimuth measurement data after outlier removal. is the threshold value corresponding to the turntable azimuth angle measurement data; when the sampling data is the turntable pitch angle measurement data, the data after outlier removal is the turntable pitch angle measurement data after outlier removal, is the threshold corresponding to the turntable pitch angle measurement data; S22. Set the :Assume that the angular acceleration of the turntable azimuth or turntable pitch angle is , then: (2) In formula (2), represents the sampling interval, y n-1 Represents the sampling data corresponding to the n-1th moment, y n-2 Represents the sampling data corresponding to the n-2th moment; S23. Combining equations (1) and (2), the turntable azimuth angle measurement data and turntable pitch angle measurement data obtained in step S1 are respectively subjected to outlier elimination using the improved limiting filter algorithm. The formulas are as follows: (3) In formula (3), Indicates The sampling data corresponding to the moment, and the other parameters are the same as formula (1).

[0014] Furthermore, the step S3 specifically comprises the following steps: S31. Determine the state model for the target angle, that is, the turntable azimuth angle or turntable pitch angle as follows: (4) In formula (4), Indicates The system state matrix at each moment is: Indicates The system state matrix at each moment is: is a natural number greater than or equal to 1, Indicates The state transfer matrix at the moment is: Indicates The control input matrix at each moment is: Indicates The amount of control at a moment, represents process noise, and its expected value is 0; S32. Determine the observation model for the target angle as: (5) In formula (5), Indicates The observed value or measured value of the system state matrix at each moment corresponds to the target angle measurement data after outliers are removed, that is, Indicates The target angle measurement data after outliers are removed at the moment, that is, The turntable azimuth measurement data after outliers are eliminated or the turntable pitch angle measurement data after outliers are eliminated; represents the state observation matrix, Indicates The measurement noise at each moment; S33. Determine the prediction equation for the target angle as: (6) The prediction equation shown in formula (6) is obtained based on the state model in step S31, where: Indicates The prior state estimate at the moment, Indicates The optimal estimate at the moment, that is, The target angle corresponding to the moment is estimated as the real data, Indicates The prior estimated covariance at time , Indicates The posterior estimated covariance at the moment, that is, The optimal estimated covariance at time, A T represents the state transfer matrix, represents the system noise covariance matrix, and the other parameters are the same as in formula (4).

[0015] S34. Based on the iterative initial value, the target angle measurement data after outliers are eliminated, and the prediction equation, iteratively calculate the optimal estimate value at each moment to obtain the target angle estimation real data corresponding to each moment; wherein the iterative initial value includes the preset initial optimal estimate value and the initial posterior estimated covariance .

[0016] Furthermore, the step S34 may include the following steps: a1. Determine the constant state transfer matrix , control input matrix and the state observation matrix , determine the system noise covariance matrix , determine the observation noise covariance , together as estimated parameters; a2. Determine the optimal estimate corresponding to the current iteration for the target angle and the posterior estimated covariance , together constitute the input value of the current iteration; among them, the optimal estimate corresponding to the first iteration and the posterior estimated covariance are the preset initial optimal estimates and the initial posterior estimated covariance ; a3. Based on the estimated parameters, the input value of the current iteration and the prediction equation, the prior state estimate and the prior estimate covariance of the target angle are calculated to obtain the following results: (7) a4. Based on the calculated prior estimated covariance and the observation noise covariance , use the Kalman gain calculation formula to calculate the Kalman gain ; The calculation formula of the Kalman gain is: (8) In formula (8), Indicates The Kalman gain at the moment, H T represents the state observation matrix; a5. Based on the calculated Kalman gain , prior state estimate and The target angle measurement data after outliers are removed at each moment , using the optimal estimate calculation formula, calculate the optimal estimate ; The optimal estimate calculation formula is: (9) In formula (9), the optimal estimate is Indicates The target angle estimation real data at each moment; a6. Based on the calculated prior estimated covariance and Kalman gain , using the optimal estimation covariance calculation formula, calculate the posterior estimation covariance ; Wherein, the optimal estimated covariance calculation formula is: (10) In formula (10), the posterior estimated covariance Indicates The optimal estimated covariance at each moment; a7. The calculated optimal estimate and the posterior estimated covariance Together they constitute the input value of a new iteration, and the iteration is repeated until the true target angle estimation data corresponding to each moment is calculated.

[0017] Furthermore, the step S34 may also include the following steps: b1. Determine the state transfer matrix , control input matrix , the state observation matrix , determine the system noise covariance matrix according to the control system accuracy , the observation noise covariance is determined according to the observation model and sensor accuracy , together as estimated parameters; b2. Determine the optimal estimate corresponding to the current iteration for the target angle and the posterior estimated covariance , together constitute the input value of the current iteration; Among them, the optimal estimate corresponding to the first iteration is and the posterior estimated covariance are the preset initial optimal estimates and the initial posterior estimated covariance ; b3. Calculate the prior state estimate of the target angle based on the estimated parameters, the input value of the current iteration and the updated prediction equation and the a priori estimated covariance ; Wherein, the updated prediction equation is: (11); b4. Based on the calculated prior estimated covariance and the observation noise covariance , the Kalman gain is calculated using the updated Kalman gain calculation formula ; The updated Kalman gain calculation formula is: (12); b5. Based on the calculated Kalman gain , prior state estimate and The target angle measurement data after outliers are removed at each moment , using the updated optimal estimate calculation formula, calculate the optimal estimate , indicating the The real data of target angle estimation at each moment; The updated optimal estimate value calculation formula is: (13); b6. Based on the calculated prior estimated covariance and Kalman gain , using the updated optimal estimation covariance calculation formula, calculate the posterior estimation covariance , indicating the The optimal estimated covariance at each moment; The updated optimal estimated covariance calculation formula is: (14) b7. The calculated optimal estimate and the posterior estimated covariance Together they constitute the input value of a new iteration, and the iteration is repeated until the true target angle estimation data corresponding to each moment is calculated.

[0018] Furthermore, the formula used in step S4 is: (15) In formula (15), Indicates The actual similarity calculated at each moment is Indicates The theoretical data of the turntable azimuth at a moment, Indicates The real data of the turntable azimuth estimation at the moment, Indicates The theoretical data of the turntable pitch angle at a certain moment, Indicates The actual data of the turntable pitch angle estimation at the moment, is a natural number.

[0019] Furthermore, in the step S5, the similarity function is fitted using the least square method according to all the moments and the corresponding actual similarities.

[0020] Furthermore, in the step S6, for each of the N consecutive new moments, steps S1 to S4 are executed to obtain the actual similarity of the moment; steps S1 to S5 are executed to obtain the estimated similarity of the moment; and then the absolute value of the difference between the actual similarity and the estimated similarity at the moment is calculated as the similarity difference value.

[0021] Furthermore, in step S7, judging whether the similarity difference values ​​of the N consecutive new moments meet the change condition includes: judging whether among the similarity difference values ​​of the N consecutive new moments, whether M similarity difference values ​​are all greater than a preset threshold, and if so, determining that the change condition is met; wherein M is a natural number greater than 0, and N / 2 <M<N。

[0022] Compared with the prior art, the present invention has the following beneficial effects: 1. The present invention tracks the changes in the azimuth and pitch angle of the turntable during the imaging process, establishes a mathematical model, evaluates the gap between the measured data and the theoretical data to obtain the actual similarity, and performs linear function fitting to predict the estimated similarity at the new moment. The difference between the estimated similarity at the new moment and the actual similarity is used to quickly discover the abnormal movement of the space target. This method can improve the timeliness of research and judgment and provide support for data application.

[0023] 2. Compared with the traditional motion detection method based on the position error of the target in the camera field of view, the present invention does not need to complete astronomical positioning, the method is simpler, and can quickly detect the target orbit movement while performing characteristic imaging of the target, with higher timeliness and realization of resource integration and utilization.

[0024] 3. The method of the present invention fills the gap in detecting and distinguishing target orbital movements using small-aperture telescopes on on-orbit satellites. It can quickly detect target orbital movements based on imaging auxiliary data, effectively reducing the risk of space target collision caused by orbital maneuvers, and can improve the safety of space target on-orbit operation. BRIEF DESCRIPTION OF THE DRAWINGS

[0025] Figure 1 A schematic flow chart of a method for distinguishing space target movement based on imaging auxiliary data provided by an embodiment of the present invention.

[0026] Figure 2 A schematic diagram of an iterative process for obtaining real data of target angle estimation according to an embodiment of the present invention.

[0027] Figure 3 A schematic diagram of another iterative process for obtaining real data of target angle estimation according to an embodiment of the present invention.

[0028] Figure 4 A process understanding diagram of a method for distinguishing space target movement based on imaging auxiliary data provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0029] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0030] See also Figure 1 , Step 1: For each moment, calculate the theoretical data of the turntable azimuth angle and the theoretical data of the turntable pitch angle of the space target according to the target orbit, and obtain the turntable azimuth angle measurement data and the turntable pitch angle measurement data of the space target during the tracking and imaging process.

[0031] The space target of the embodiment of the present invention may be various types of spacecraft in space, and is not specifically limited here.

[0032] Specifically, step 1 includes: using the orbit of the space target and the azimuth and pitch measurement values ​​of the turntable to calculate the theoretical data of the turntable azimuth angle and the theoretical data of the turntable pitch angle of the space target at each moment, and obtaining the turntable azimuth angle measurement data of the space target at each moment output in real time by the turntable azimuth encoder, and the turntable pitch angle measurement data of the space target at each moment output in real time by the turntable pitch encoder.

[0033] It should be noted that the turntable azimuth encoder and turntable pitch encoder are devices in the process of space target tracking and imaging.

[0034] Through step 1, for each moment, theoretical data of the turntable azimuth angle, theoretical data of the turntable pitch angle, measurement data of the turntable azimuth angle and measurement data of the turntable pitch angle at that moment can be obtained, which belong to the imaging auxiliary data in the embodiment of the present invention.

[0035] Step 2: Use the improved limiting filter algorithm to remove outliers from the obtained turntable azimuth angle measurement data and turntable pitch angle measurement data.

[0036] Since there may be outliers in the obtained turntable azimuth angle measurement data and turntable pitch angle measurement data, which will have a great impact on the abnormal movement judgment of the space target, the embodiment of the present invention uses an improved limiting filtering algorithm to eliminate outliers for the obtained turntable azimuth angle measurement data and turntable pitch angle measurement data.

[0037] The embodiment of the present invention first proposes the following initial formula for outlier removal based on the limiting filtering algorithm: (1) In formula (1), represents the data after outliers are removed. Indicates The sampling data corresponding to the moment, Indicates The sampling data corresponding to the moment, is a natural number greater than or equal to 1.

[0038] Indicates the threshold value corresponding to the sampling data. When the sampling data is the turntable azimuth measurement data, the data after outlier removal is the turntable azimuth measurement data after outlier removal. is the threshold value corresponding to the turntable azimuth angle measurement data; when the sampling data is the turntable pitch angle measurement data, the data after outlier removal is the turntable pitch angle measurement data after outlier removal, is the threshold corresponding to the turntable pitch angle measurement data.

[0039] For the turntable azimuth measurement data and turntable pitch angle measurement data, the corresponding threshold Different, for example, for the turntable azimuth measurement data and the turntable pitch angle measurement data, the corresponding thresholds They may be the maximum angular velocity of the turntable azimuth angle and the maximum angular velocity of the turntable pitch angle, respectively.

[0040] By setting , for each current moment The corresponding turntable azimuth angle measurement data and turntable pitch angle measurement data are respectively eliminated for outliers to obtain the corresponding .

[0041] To set , assuming that the angular acceleration of the turntable azimuth or turntable pitch angle is , then: (2) In formula (2), represents the sampling interval, y n-1 Represents the sampling data corresponding to the n-1th moment, y n-2 Represents the sampling data corresponding to the n-2th moment.

[0042] Then, the obtained turntable azimuth angle measurement data and turntable pitch angle measurement data are respectively subjected to outlier elimination using the improved limiting filter algorithm, and the formula used can be: (3) In formula (3), Indicates The sampling data corresponding to the moment, and the other parameters are the same as formula (1).

[0043] The traditional method does not pre-process the turntable azimuth measurement data and the turntable pitch angle measurement data which are the original data. The wild values ​​in these original data may lead to the misjudgment of the movement of the space target. At the same time, the error between the measured value and the true value may make it impossible to complete the judgment. The embodiment of the present invention proposes an improved limiting filtering algorithm, which pre-processes the original data to remove wild values, can improve the data processing efficiency and accuracy, lay a solid foundation for subsequent data processing, can effectively avoid the misjudgment caused by wild values, and can improve the accuracy of the space target movement judgment.

[0044] Step 3: Perform filtering estimation on the turntable azimuth angle measurement data and the turntable pitch angle measurement data after outliers are eliminated, respectively, to obtain the turntable azimuth angle estimation true data and the turntable pitch angle estimation true data corresponding to each moment.

[0045] In an optional implementation manner, step 3 may include the following steps 31 to 34.

[0046] In the embodiment of the present invention, steps 31 to 34 are respectively performed for the turntable azimuth angle and the turntable pitch angle to obtain the turntable azimuth angle estimation real data and the turntable pitch angle estimation real data corresponding to each moment.

[0047] In step 31 to step 34, the description is made in terms of a target angle, where the target angle is a turntable azimuth angle or a turntable pitch angle.

[0048] Step 31: Determine the state model for the target angle, that is, the turntable azimuth angle or the turntable pitch angle as follows: (4) In formula (4), Indicates The system state matrix at each moment is: Indicates The system state matrix at each moment is: is a natural number greater than or equal to 1, Indicates The state transfer matrix at the moment is: Indicates The control input matrix at each moment is: Indicates The amount of control at a moment, represents process noise, and its expected value is 0.

[0049] Step 32: Determine the observation model for the target angle as: (5) In formula (5), Indicates The observed value or measured value of the system state matrix at each moment corresponds to the target angle measurement data after outliers are removed, that is, Indicates The target angle measurement data after outliers are removed at the moment, that is, The turntable azimuth measurement data after outliers are eliminated or the turntable pitch angle measurement data after outliers are eliminated; represents the state observation matrix, Indicates The measurement noise at that moment.

[0050] Step 33: Determine the prediction equation for the target angle as: (6) The prediction equation shown in formula (6) is obtained based on the state model in step S31, where: Indicates The prior state estimate at the moment, Indicates The optimal estimate at the moment, that is, The target angle corresponding to the moment is estimated as the real data, Indicates The prior estimated covariance at time , Indicates The posterior estimated covariance at the moment, that is, The optimal estimated covariance at time, A T represents the state transfer matrix, represents the system noise covariance matrix, and the other parameters are the same as in formula (4).

[0051] Step 34: Based on the iterative initial value, the target angle measurement data after outliers are eliminated, and the prediction equation, iteratively calculate the optimal estimate value at each moment to obtain the target angle estimated real data corresponding to each moment; wherein the iterative initial value includes the preset initial optimal estimate value and the initial posterior estimated covariance .

[0052] For details, see Figure 2 , step S34 may include the following steps: Step a1: Determine the constant state transfer matrix , control input matrix and the state observation matrix , determine the system noise covariance matrix , determine the observation noise covariance , together as the estimated parameters.

[0053] Step a2: Determine the optimal estimated value corresponding to the current iteration for the target angle and the posterior estimated covariance , together constitute the input value of the current iteration; among them, the optimal estimate corresponding to the first iteration and the posterior estimated covariance are the preset initial optimal estimates and the initial posterior estimated covariance .

[0054] Step a3: Based on the estimated parameters, the input value of the current iteration and the prediction equation, the prior state estimation value and the prior estimation covariance of the target angle are calculated to obtain the following results: (7).

[0055] Step a4: Estimate the covariance based on the calculated prior and the observation noise covariance , use the Kalman gain calculation formula to calculate the Kalman gain ; The calculation formula of the Kalman gain is: (8) In formula (8), Indicates The Kalman gain at the moment, H T represents the state observation matrix.

[0056] Step a5: Based on the calculated Kalman gain , prior state estimate and The target angle measurement data after outliers are removed at each moment , using the optimal estimate calculation formula, calculate the optimal estimate ; The optimal estimate calculation formula is: (9) In formula (9), the optimal estimate is Indicates The target angle estimation real data at each moment.

[0057] Step a6: Estimate the covariance based on the calculated prior and Kalman gain , using the optimal estimation covariance calculation formula, calculate the posterior estimation covariance ; Wherein, the optimal estimated covariance calculation formula is: (10) In formula (10), the posterior estimated covariance Indicates The optimal estimated covariance at each moment.

[0058] It is understandable that for At this moment, through the above steps a1~a6, the optimal estimate can be obtained As the The target angle estimation real data at each moment , and obtain the posterior estimated covariance , which is equivalent to completing the The actual data estimate at that moment.

[0059] Step a7: The calculated optimal estimate and the posterior estimated covariance Together they constitute the input value of a new iteration, and the iteration is repeated until the true target angle estimation data corresponding to each moment is calculated.

[0060] The best estimate calculated and the posterior estimated covariance The input values ​​of the new iteration are executed from step a2 to step a6 to obtain the real data of the target angle estimation at the next moment. , and obtain the posterior estimated covariance , which is equivalent to completing the The actual data estimate at that moment.

[0061] Then, through continuous iteration, the true data of the turntable azimuth angle estimation and the true data of the turntable pitch angle estimation at each moment corresponding to step 1 can be obtained.

[0062] The embodiment of the present invention selects a suitable initial optimal estimate value and the initial posterior estimated covariance As the initial value of the iteration, the initial value of the iteration and the estimated parameter are substituted into the filter equation formed by equations (7) to (10) for iteration, so that the measured value is closer to the true value. It can be understood that the above filter equation is implemented based on Kalman filtering.

[0063] In order to be more in line with the actual situation, the embodiment of the present invention specifically determines the value of the estimation parameter.

[0064] Among them, the ideal physical model of the turntable azimuth angle and turntable pitch angle is as follows: (11) In formula (11), Indicates The turntable azimuth or turntable pitch angle at a moment, Indicates The angular velocity of the turntable azimuth or turntable pitch angle at a certain moment, represents the angular acceleration, Indicates the sampling period.

[0065] Converting equation (11) into matrix form, we have: (12) make , then: (13) In the actual system, the satellite remote sensing data only contains angle information, not angular velocity and angular acceleration information, so let the angular velocity , angular acceleration , so formula (12) can be simplified as: (14) In formula (14), Represents a selected value, the value range is .

[0066] Therefore, it can be concluded that in the actual system: the state transfer matrix , control input matrix , the system noise covariance matrix The specific value is determined according to the accuracy of the control system. The parameters obtained by observation are still the turntable azimuth angle and turntable pitch angle. Therefore, the observation model corresponding to equation (5) can be changed to: (15) in, Represents a selected value, the value range is

[0067] Therefore, it can be concluded that in the actual system: the state observation matrix , observation noise covariance The specific value is determined according to the observation model and the sensor accuracy. Therefore, compared with step a1 to step a7, the filter equation can be simplified for the actual system.

[0068] For details, see Figure 3 Step 34 may include the following steps b1 to b7: Step b1: Determine the state transfer matrix , control input matrix , the state observation matrix , determine the system noise covariance matrix according to the control system accuracy , the observation noise covariance is determined according to the observation model and sensor accuracy , together as the estimated parameters.

[0069] Step b2: Determine the optimal estimated value corresponding to the current iteration for the target angle and the posterior estimated covariance , together constitute the input value of the current iteration; Among them, the optimal estimate corresponding to the first iteration is and the posterior estimated covariance are the preset initial optimal estimates and the initial posterior estimated covariance .

[0070] Step b3: Calculate the prior state estimate of the target angle based on the estimated parameters, the input value of the current iteration and the updated prediction equation and the a priori estimated covariance ; Wherein, the updated prediction equation is: (16).

[0071] Step b4: Estimate the covariance based on the calculated prior and the observation noise covariance , the Kalman gain is calculated using the updated Kalman gain calculation formula ; The updated Kalman gain calculation formula is: (17).

[0072] Step b5: Based on the calculated Kalman gain , prior state estimate and The target angle measurement data after outliers are removed at each moment , using the updated optimal estimate calculation formula, calculate the optimal estimate , indicating the The real data of target angle estimation at each moment; The updated optimal estimate value calculation formula is: (18).

[0073] Step b6: Estimate the covariance based on the calculated prior and Kalman gain , using the updated optimal estimation covariance calculation formula, calculate the posterior estimation covariance , indicating the The optimal estimated covariance at each moment; The updated optimal estimated covariance calculation formula is: (19).

[0074] It is understandable that for At this moment, through the above steps b1 to b6, the optimal estimate can be obtained As the The target angle estimation real data at each moment , and obtain the posterior estimated covariance Indicates The optimal estimated covariance at the moment completes the The actual data estimate at that moment.

[0075] Step b7: The calculated optimal estimate and the posterior estimated covariance Together they constitute the input value of a new iteration, and the iteration is repeated until the true target angle estimation data corresponding to each moment is calculated.

[0076] The best estimate calculated and the posterior estimated covariance The input values ​​of the new iteration are executed from step b2 to step b6 to obtain the target angle estimation real data at the next moment. , and obtain the posterior estimated covariance , which is equivalent to completing the The actual data estimate at that moment.

[0077] Then, through continuous iteration, the true data of the turntable azimuth angle estimation and the true data of the turntable pitch angle estimation at each moment corresponding to step 1 can be obtained.

[0078] In the embodiment of the present invention, the state transition matrix is ​​determined , control input matrix , the state observation matrix , determine the system noise covariance matrix according to the control system accuracy , the observation noise covariance is determined according to the observation model and sensor accuracy , it is possible to obtain estimated parameters that conform to the actual system, and to make the true data of the turntable azimuth angle estimation and the turntable pitch angle estimation closer to the true data, thereby improving the accuracy of subsequent space target movement discrimination.

[0079] The embodiment of the present invention can input measurement values, namely, turntable azimuth angle measurement data and turntable pitch angle measurement data, and can effectively estimate the true values ​​of the turntable azimuth angle and pitch angle based on the Kalman filter algorithm, realize the correction of the measurement values, and reduce the relative error between the measurement value and the true value.

[0080] Step 4: Using the theoretical data of the turntable azimuth angle, the estimated real data of the turntable azimuth angle, the theoretical data of the turntable pitch angle and the estimated real data of the turntable pitch angle at each moment, calculate the actual similarity between the estimated real data and the theoretical data at the corresponding moment.

[0081] The formula used in step 4 is: (20) In formula (20), Indicates The actual similarity calculated at each moment is Indicates The theoretical data of the turntable azimuth at a moment, Indicates The real data of the turntable azimuth estimation at the moment, Indicates The theoretical data of the turntable pitch angle at a certain moment, Indicates The actual data of the turntable pitch angle estimation at the moment, is a natural number.

[0082] It can be understood that formula (20) is implemented based on the Euclidean distance calculation formula. , the corresponding actual similarity can be calculated Therefore, multiple sets of data can be obtained ,in, Indicates time .

[0083] Step 5: Fit the similarity function based on the calculated actual similarities at all times; In an optional implementation manner, step 5 may include: According to all the moments and their corresponding actual similarities, the similarity function is fitted using the least squares method to predict the estimated similarity of the new moment.

[0084] Specifically, according to the multiple sets of data obtained in step 4 The least squares method can be used to fit a linear equation to represent the similarity function, which can be specifically: (twenty one) In formula (21), Represents the estimated similarity, that is, by inputting a moment into the similarity function , the predicted The similarity between the estimated real data and the theoretical data at the moment, and The fitting coefficients can be determined by the least squares method.

[0085] The embodiment of the present invention uses the actual similarities of multiple historical moments to fit a similarity function, and can predict the corresponding estimated similarity for each new moment, so as to subsequently use the difference between the estimated similarities of multiple new moments and the actual similarities obtained by actual calculation to determine the changes in the spatial target.

[0086] Step 6: for each of the N consecutive new moments, obtain the actual similarity of the moment, use the similarity function to obtain the estimated similarity of the moment, and calculate the absolute value of the difference between the actual similarity and the estimated similarity at the moment as the similarity difference value; Wherein, N is a natural number greater than 0, which can be selected as needed, for example, 8, 10, etc.

[0087] It should be noted that the N consecutive new moments in step 6 refer to moments for which actual similarity has not been calculated before, nor have they been used to fit the similarity function, and no spatial target movement discrimination has been performed. For the first execution of step 6, the N consecutive new moments refer to multiple sets of data for fitting the similarity function. The next N consecutive new moments.

[0088] For each moment in the N consecutive new moments, the actual similarity of the moment is obtained, which can be performed according to steps 1 to 4, and will not be repeated here.

[0089] For each of the N consecutive new moments, the moment can be substituted into the similarity function shown in formula (21) to obtain the estimated similarity of the moment.

[0090] Then, the absolute value of the difference between the actual similarity and the estimated similarity at that moment is calculated as the similarity difference value.

[0091] Step 7: Determine whether the similarity difference values ​​of the N consecutive new moments meet the change condition. If so, determine that the space target has changed; if not, return to step 6.

[0092] In an optional implementation, the determining whether the similarity difference values ​​of the N consecutive new moments satisfy the change condition comprises: determining whether among the similarity difference values ​​of the N consecutive new moments, whether M similarity difference values ​​are all greater than a preset threshold, and if so, determining that the change condition is satisfied; wherein M is a natural number greater than 0, and N / 2 <M<N。

[0093] It is understandable that satisfying the change condition requires that more than half of the similarity difference values ​​of the N consecutive new moments are greater than a preset threshold, which can avoid misjudgment.

[0094] In the specific implementation process, each of the N consecutive new moments can be judged in turn. If the similarity difference value of the current moment is greater than the preset threshold, the counter is increased by one; when the counter count value reaches M, it is determined that the abnormal condition is met, and it is determined that the spatial target has an abnormal movement; otherwise, return to step 6 to continue monitoring and abnormal movement judgment for each of the N consecutive new moments that appear again.

[0095] As a preferred implementation, N can be 10 and M can be 6, that is, the similarity difference value calculation results are judged as 6 out of 10, that is, if 6 out of 10 consecutive similarity difference values ​​are greater than the preset threshold, it is determined that the space target has moved, otherwise return to step 6 to continue monitoring and judging.

[0096] The implementation process of the method for distinguishing space target movement based on imaging auxiliary data provided by the embodiment of the present invention is shown in Figure 4 For details, please refer to the relevant steps in the previous article and will not be repeated here.

[0097] In the embodiment of the present invention, a similarity function is fitted through a large amount of data of a space target that has not changed its orbit, and a similarity function is used to perform state estimation to evaluate whether the space target has changed its orbit. Specifically, a corresponding similarity difference value is obtained by comparing the actual similarity and the estimated similarity at each moment, and compared with a preset threshold. When the number of moments when the similarity difference value is greater than the preset threshold reaches a certain number, it is determined that the space target has moved. This method does not require astronomical positioning, and the method is simpler. It can quickly detect the target orbit movement while performing characteristic imaging on the target, and has higher timeliness.

[0098] In the method for distinguishing the abnormal movement of a space target based on imaging auxiliary data provided by an embodiment of the present invention, for the azimuth and pitch angle of the turntable of the space target, firstly, theoretical data and measurement data at each moment are obtained; then, the measured data are subjected to outlier elimination using an improved limiting filter algorithm to improve data processing efficiency and accuracy, and avoid misjudgment caused by outliers; next, the measured data after the outliers are eliminated are subjected to filtering estimation to obtain estimated real data corresponding to each moment to correct the measured data and reduce relative errors; then, the actual similarity between the estimated real data and the theoretical data at each moment is calculated; and a similarity function is fitted according to the calculated actual similarities of all moments; then, for each moment in N consecutive new moments, the actual similarity of the moment is obtained, the estimated similarity of the moment is obtained using the similarity function, the absolute value of the difference between the actual similarity and the estimated similarity at the moment is calculated as the similarity difference value, and whether the space target has moved is judged by judging whether the similarity difference values ​​of the N consecutive new moments meet the abnormal movement condition.

[0099] The embodiment of the present invention tracks the changes in the azimuth and pitch angle of the turntable during the imaging process, establishes a mathematical model, evaluates the gap between the measured data and the theoretical data to obtain the actual similarity, and performs linear function fitting to predict the estimated similarity at the new moment, and uses the difference between the estimated similarity at the new moment and the actual similarity to quickly discover the movement of the space target. This method can improve the timeliness of research and judgment and provide support for data application. In addition, compared with the traditional movement detection method based on the position error of the target in the camera field of view, the present invention does not need to complete astronomical positioning, the method is simpler, and can quickly detect the movement of the target orbit while performing characteristic imaging of the target, with higher timeliness, and realizes the integration and application of resources.

[0100] The method of the present invention fills the blank of detecting and distinguishing target orbital movements using small-aperture telescopes on on-orbit satellites. It can quickly detect target orbital movements based on imaging auxiliary data, effectively reduce the risk of space target collision caused by orbital maneuvers, and improve the safety of space target on-orbit operation.

[0101] Although this specification is described according to implementation modes, not every implementation mode includes only one independent technical solution. This description of the specification is only for the sake of clarity. Those skilled in the art should regard the specification as a whole. The technical solutions in each embodiment may also be appropriately combined to form other implementation modes that can be understood by those skilled in the art.

[0102] Therefore, the above description is only a preferred embodiment of the present application and is not intended to limit the scope of implementation of the present application; that is, all equivalent changes made according to the scope of the claims of the present application are within the protection scope of the claims of the present application.

Claims

1. A method for distinguishing space target movement based on imaging auxiliary data, characterized in that: The specific steps include: S1. At each moment, calculate the theoretical data of the turntable azimuth angle and the theoretical data of the turntable pitch angle of the space target according to the target orbit, and obtain the turntable azimuth angle measurement data and the turntable pitch angle measurement data of the space target during the tracking and imaging process; S2. The obtained turntable azimuth angle measurement data and turntable pitch angle measurement data are respectively subjected to outlier elimination using an improved limiting filter algorithm; S3. Perform filtering estimation on the turntable azimuth angle measurement data and the turntable pitch angle measurement data after outliers are eliminated, respectively, to obtain the turntable azimuth angle estimation true data and the turntable pitch angle estimation true data corresponding to each moment; S4. Using the obtained theoretical data of the turntable azimuth angle, the estimated real data of the turntable azimuth angle, the theoretical data of the turntable pitch angle and the estimated real data of the turntable pitch angle at each time, calculate the actual similarity between the estimated real data and the theoretical data at the corresponding time; S5. According to the calculated actual similarities of all moments, a similarity function is fitted to calculate the estimated similarity of the moment; S6. for each of the N consecutive new moments, obtain the actual similarity of the moment, use the similarity function to obtain the estimated similarity of the moment, and calculate the absolute value of the difference between the actual similarity and the estimated similarity at the moment as the similarity difference value, where N is a natural number greater than 0; S7. Determine whether the similarity difference values ​​of the N consecutive new moments meet the change condition, if so, determine that the space target has changed; if not, return to step S6.

2. The method for distinguishing space target movement based on imaging auxiliary data according to claim 1, characterized in that: The step S1 specifically includes: using the orbit of the space target and the azimuth and pitch measurement values ​​of the turntable to calculate the theoretical data of the turntable azimuth angle and the theoretical data of the turntable pitch angle of the space target at each moment, and obtaining the turntable azimuth angle measurement data of the space target at each moment output in real time by the turntable azimuth encoder, and the turntable pitch angle measurement data of the space target at each moment output in real time by the turntable pitch encoder.

3. The method for distinguishing space target movement based on imaging auxiliary data according to claim 1, characterized in that: The step S2 specifically includes the following steps: S21. Based on the limiting filter algorithm, the following initial formula for outlier removal is proposed: （1） In formula (1), represents the data after outliers are removed. Indicates The sampling data corresponding to the moment, Indicates The sampling data corresponding to the moment, is a natural number greater than or equal to 1; Indicates the threshold value corresponding to the sampling data. When the sampling data is the turntable azimuth measurement data, the data after outlier removal is the turntable azimuth measurement data after outlier removal. is the threshold value corresponding to the turntable azimuth angle measurement data; when the sampling data is the turntable pitch angle measurement data, the data after outlier removal is the turntable pitch angle measurement data after outlier removal, is the threshold corresponding to the turntable pitch angle measurement data; S22. Set the :Assume that the angular acceleration of the turntable azimuth or turntable pitch angle is , then: （2） In formula (2), represents the sampling interval, y n-1 Represents the sampling data corresponding to the n-1th moment, y n-2 Represents the sampling data corresponding to the n-2th moment; S23. Combining equations (1) and (2), the turntable azimuth angle measurement data and turntable pitch angle measurement data obtained in step S1 are respectively subjected to outlier elimination using the improved limiting filter algorithm. The formulas are as follows: （3） In formula (3), Indicates The sampling data corresponding to the moment, and the other parameters are the same as formula (1).

4. The method for distinguishing space target movement based on imaging auxiliary data according to claim 1, characterized in that: The step S3 specifically comprises the following steps: S31. Determine the state model for the target angle, that is, the turntable azimuth angle or turntable pitch angle as follows: （4） In formula (4), Indicates The system state matrix at each moment is: Indicates The system state matrix at each moment is: is a natural number greater than or equal to 1, Indicates The state transfer matrix at the moment is: Indicates The control input matrix at each moment is: Indicates The amount of control at a moment, represents process noise, and its expected value is 0; S32. Determine the observation model for the target angle as: （5） In formula (5), Indicates The observed value or measured value of the system state matrix at each moment corresponds to the target angle measurement data after outliers are removed, that is, Indicates The target angle measurement data after outliers are removed at the moment, that is, The turntable azimuth measurement data after outliers are eliminated or the turntable pitch angle measurement data after outliers are eliminated; represents the state observation matrix, Indicates The measurement noise at each moment; S33. Determine the prediction equation for the target angle as: （6） The prediction equation shown in formula (6) is obtained based on the state model in step S31, where: Indicates The prior state estimate at the moment, Indicates The optimal estimate at the moment, that is, The target angle estimation real data corresponding to the moment, Indicates The prior estimated covariance at time , Indicates The posterior estimated covariance at the moment, that is, The optimal estimated covariance at time, A T represents the state transfer matrix, represents the system noise covariance matrix, and the other parameters are the same as in formula (4); S34. Based on the iterative initial value, the target angle measurement data after outliers are eliminated, and the prediction equation, iteratively calculate the optimal estimate value at each moment to obtain the target angle estimation real data corresponding to each moment; wherein the iterative initial value includes the preset initial optimal estimate value and the initial posterior estimated covariance .

5. The method for distinguishing space target movement based on imaging auxiliary data according to claim 4, characterized in that: The step S34 may include the following steps: a1. Determine the constant state transfer matrix , control input matrix and the state observation matrix , determine the system noise covariance matrix , determine the observation noise covariance , together as estimated parameters; a2. Determine the optimal estimate corresponding to the current iteration for the target angle and the posterior estimated covariance , together constitute the input value of the current iteration; among them, the optimal estimate corresponding to the first iteration and the posterior estimated covariance are the preset initial optimal estimates and the initial posterior estimated covariance ; a3. Based on the estimated parameters, the input value of the current iteration and the prediction equation, the prior state estimate and the prior estimate covariance of the target angle are calculated to obtain the following results: （7） a4. Based on the calculated prior estimated covariance and the observation noise covariance , use the Kalman gain calculation formula to calculate the Kalman gain ; The calculation formula of the Kalman gain is: （8） In formula (8), Indicates The Kalman gain at the moment, H T represents the state observation matrix; a5. Based on the calculated Kalman gain , prior state estimate and The target angle measurement data after outliers are removed at each moment , using the optimal estimate calculation formula, calculate the optimal estimate ; The optimal estimate calculation formula is: （9） In formula (9), the optimal estimate is Indicates The target angle estimation real data at each moment; a6. Based on the calculated prior estimated covariance and Kalman gain , using the optimal estimation covariance calculation formula, calculate the posterior estimation covariance ; Wherein, the optimal estimated covariance calculation formula is: （10） In formula (10), the posterior estimated covariance Indicates The optimal estimated covariance at each moment; a7. The calculated optimal estimate and the posterior estimated covariance Together they constitute the input value of a new iteration, and the iteration is repeated until the true target angle estimation data corresponding to each moment is calculated.

6. The method for distinguishing space target movement based on imaging auxiliary data according to claim 4, characterized in that: The step S34 may also include the following steps: b1. Determine the state transfer matrix , control input matrix , the state observation matrix , determine the system noise covariance matrix according to the control system accuracy , the observation noise covariance is determined according to the observation model and sensor accuracy , together as estimated parameters; b2. Determine the optimal estimate corresponding to the current iteration for the target angle and the posterior estimated covariance , together constitute the input value of the current iteration; Among them, the optimal estimate corresponding to the first iteration is and the posterior estimated covariance are the preset initial optimal estimates and the initial posterior estimated covariance ; b3. Calculate the prior state estimate of the target angle based on the estimated parameters, the input value of the current iteration and the updated prediction equation and the a priori estimated covariance ; Wherein, the updated prediction equation is: （11）； b4. Based on the calculated prior estimated covariance and the observation noise covariance , the Kalman gain is calculated using the updated Kalman gain calculation formula ; The updated Kalman gain calculation formula is: （12）； b5. Based on the calculated Kalman gain , prior state estimate and The target angle measurement data after outliers are removed at each moment , using the updated optimal estimate calculation formula, calculate the optimal estimate , indicating the The real data of target angle estimation at each moment; The updated optimal estimate value calculation formula is: （13）； b6. Based on the calculated prior estimated covariance and Kalman gain , using the updated optimal estimate covariance calculation formula, calculate the posterior estimate covariance , indicating the The optimal estimated covariance at each moment; The updated optimal estimated covariance calculation formula is: （14）； b7. The calculated optimal estimate and the posterior estimated covariance Together they constitute the input value of a new iteration, and the iteration is repeated until the true target angle estimation data corresponding to each moment is calculated.

7. The method for distinguishing space target movement based on imaging auxiliary data according to claim 1, characterized in that: The formula used in step S4 is: （15） In formula (15), Indicates The actual similarity calculated at each moment is Indicates The theoretical data of the turntable azimuth at a moment, Indicates The real data of the turntable azimuth estimation at the moment, Indicates The theoretical data of the turntable pitch angle at a certain moment, Indicates The actual data of the turntable pitch angle estimation at the moment, is a natural number.

8. The method for distinguishing space target movement based on imaging auxiliary data according to claim 1, characterized in that: In the step S5, the similarity function is fitted using the least square method according to all the moments and the corresponding actual similarities.

9. The method for distinguishing space target movement based on imaging auxiliary data according to claim 1, characterized in that: In step S6, for each of the N consecutive new moments, steps S1 to S4 are executed to obtain the actual similarity of the moment; steps S1 to S5 are executed to obtain the estimated similarity of the moment; and then the absolute value of the difference between the actual similarity and the estimated similarity at the moment is calculated as the similarity difference value.

10. The method for distinguishing space target movement based on imaging auxiliary data according to claim 1, characterized in that: In the step S7, judging whether the similarity difference values ​​of the N consecutive new moments meet the change condition includes: judging whether among the similarity difference values ​​of the N consecutive new moments, whether M similarity difference values ​​are all greater than a preset threshold value, and if so, determining that the change condition is met; wherein M is a natural number greater than 0, and N / 2 <M<N。

Images