A method for identifying spatial target abnormality based on imaging auxiliary data

Abstract

The application discloses a kind of space target abnormal change discrimination methods based on imaging auxiliary data, first acquire the azimuth and pitch angle of space target turntable theory and measured data, using improved amplitude limiting filtering algorithm to eliminate outliers;Filter estimation is carried out to the data after eliminating outliers, to obtain estimated real data to correct measured data;The actual similarity between the estimated real data and the theoretical data is calculated, and the similarity function is fitted;For continuous new time, the absolute value of the difference between the actual similarity and the estimated similarity obtained by the similarity function is taken as the similarity difference value.The application does not need astronomical positioning, can detect target orbit abnormal change while imaging target characteristics, improve the timeliness of research and judgment and the safety of space target on-orbit operation, effectively reduce the risk of space target collision, fill the blank of detecting target orbit abnormal change by using on-orbit satellite small aperture telescope.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of space situation awareness and space target identification, and particularly relates to a space target abnormal behavior identification method based on imaging auxiliary data. BACKGROUND

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

[0003] In the low orbit, with the deployment of a series of giant constellations, the space will be severely squeezed. Such a large number and densely distributed satellites will bring huge collision risks to the launch, normal operation and disposal of spacecraft. Starlink satellites will continue to lift the orbit for about 1 month under the action of small thrust in the initial launch into orbit, and there will also be a phenomenon of thrust value jump before and after the start and stop of thrust. At the same time, in order to avoid the collision risk within the constellation, the Starlink satellite is equipped with an autonomous collision avoidance system, which can receive space target catalog data and perform orbit maneuver when needed. Starlink satellites will often make small-scale orbit correction and adjustment to avoid collision, but Space Exploration Company does not completely release the orbit maneuver of Starlink satellites, which will increase the risk of operation of other spacecraft by several times. In the high orbit, with the development of satellite platform from large integration to small decentralization, high-orbit satellites more frequently perform orbit maneuver to carry out satellite life extension, space test and other tasks.

[0004] Currently, spacecraft in orbit maneuver frequently, especially when there is a collision warning event between domestic and foreign spacecraft, both sides often do not communicate and maneuver for "back-to-back" disposal, which has the risk of collision due to maneuver. Observation, monitoring and identification of these orbit abnormal space targets have always been the focus of research. Timely discovery and identification of the orbit abnormal behavior of such targets can effectively support space collision warning. The small-aperture telescope carried by China's on-orbit spacecraft can track the target for a long time and obtain a large amount of observation data, but due to the influence of equipment capacity, the detection capability of such small-aperture telescope is insufficient, resulting in limited number of stars in the field of view, which often cannot complete astronomical positioning measurement in the target tracking process. At the same time, the angle measurement accuracy of such telescope is poor, which is difficult to position accurately, so the traditional abnormal behavior detection method cannot be used for abnormal behavior identification, and the data is difficult to fully utilize. Therefore, an effective space target abnormal behavior identification method is urgently needed.

[0005] For the detection of the abnormal motion of the space target, the principle of the traditional abnormal motion detection is: the angle position of the target in the camera coordinate system is compared with the actual angle position, and when the detection threshold is exceeded, it is considered that the abnormal motion occurs. According to the predicted position of the target and the coordinate system information of the monitoring satellite / camera, the visibility of the target relative to the camera and the position of the target in the camera coordinate system if the target is visible can be calculated The difference between the actual position of the target in the camera coordinate system and the predicted position of the target in the camera coordinate system is denoted as Whether the target appears in the field of view of the camera after the maneuver, in the case that the target is visible, is related to 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 is less than the field of view angle and the direction of the position difference is ideal, the target will appear in the field of view after the maneuver; if the predicted position is at the center of the field of view, when is less than the field of view angle , the target will appear in the field of view after the maneuver.

[0006] The traditional abnormal motion detection adopts a method based on a characteristic quantity , that is, the measured value and the theoretical predicted value are converted to calculate the characteristic quantity according to the following formula, and when the set abnormal motion threshold is exceeded, it is determined that the abnormal motion occurs, therefore, the main influencing factors of the traditional abnormal motion detection identification include the observation error of the characteristic quantity and the theoretical calculation error of the characteristic quantity, specifically, the error is calculated by using the measured value and the theoretical value of the astronomical positioning.

[0007]

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

[0009] Although the method has the advantage of high measurement accuracy, the method is complex and has poor timeliness, and it is required to first complete the astronomical positioning, the measurement angle accuracy is high, and it is not suitable for imaging equipment. SUMMARY

[0010] The purpose of the present application is to overcome the defects and deficiencies of the prior art, and to provide a space target abnormal motion identification method based on imaging auxiliary data.

[0011] In order to achieve the above purpose, the present application provides the following technical scheme:

[0012] The application discloses a space target abnormality distinguishing method based on imaging auxiliary data.

[0013] S1. According to the target orbit, the azimuth angle and the elevation angle of the rotating platform of the space target are calculated, and the azimuth angle and the elevation angle of the rotating platform of the space target in the tracking imaging process are obtained;

[0014] S2. The obtained azimuth angle and the elevation angle of the rotating platform are subjected to wild value elimination by using an improved amplitude limiting filtering algorithm.

[0015] S3. The azimuth angle and the elevation angle of the rotating platform after wild value elimination are subjected to filtering estimation, and the estimated real data of the azimuth angle and the elevation angle of the rotating platform corresponding to each time are obtained.

[0016] S4. The actual similarity between the estimated real data and the theoretical data of each time is calculated by using the obtained azimuth angle, the estimated real data of the azimuth angle, the elevation angle and the estimated real data of the elevation angle.

[0017] S5. The similarity function is fitted according to the calculated actual similarity of all times, and the estimated similarity of the time is calculated.

[0018] S6. The actual similarity of each time in N continuous new times is obtained, the estimated similarity of each time is obtained by using the similarity function, and the absolute value of the difference between the actual similarity and the estimated similarity of each time is calculated as the similarity difference value, wherein N is a natural number greater than 0.

[0019] S7. Whether the similarity difference values of the N continuous new times meet the abnormality condition is judged, if yes, it is determined that the space target has an abnormality, and if not, the step S6 is returned.

[0020] Further, the step S1 specifically comprises the following steps: the azimuth angle and the elevation angle of the rotating platform of the space target at each time are calculated by using the orbit of the space target and the azimuth and elevation measurement values of the rotating platform, the azimuth angle measurement data of the space target at each time output by the azimuth encoder of the rotating platform in real time is obtained, and the azimuth angle measurement data of the space target at each time output by the elevation encoder of the rotating platform in real time is obtained.

[0021] Further, the step S2 specifically comprises the following steps:

[0022] S21. The following initial formula for wild value elimination is proposed on the basis of the amplitude limiting filtering algorithm:

[0023] (1)

[0024] In formula (1), denotes the data after wild value elimination, denotes the sampling data corresponding to the time point, denotes the sampling data corresponding to the time point, is a natural number greater than or equal to 1;

[0025] denotes the threshold value corresponding to the sampling data, when the sampling data is the turntable azimuth measurement data, the data after wild value elimination is the wild value eliminated turntable azimuth measurement data, is the threshold value corresponding to the turntable azimuth measurement data; when the sampling data is the turntable pitch angle measurement data, the data after wild value elimination is the wild value eliminated turntable pitch angle measurement data, is the threshold value corresponding to the turntable pitch angle measurement data;

[0026] S22. Set in formula (1) : assuming that the rotation angle acceleration of the turntable azimuth or the turntable pitch angle is , then:

[0027] (2)

[0028] In formula (2), denotes the sampling interval, y n-1 denotes the sampling data corresponding to the n-1 time point, y n-2 denotes the sampling data corresponding to the n-2 time point;

[0029] S23. In combination with formula (1) and (2), the turntable azimuth measurement data and the turntable pitch angle measurement data obtained by step S1 are respectively subjected to wild value elimination by using the improved amplitude limiting filtering algorithm, and the formula is as follows:

[0030] (3)

[0031] In formula (3), denotes the sampling data corresponding to the time point, and the remaining parameters are the same as formula (1).

[0032] Further, the step S3 specifically includes the following steps:

[0033] S31. Determine that the state model for the target angle, i.e. the turntable azimuth or the turntable pitch angle, is:

[0034] (4)

[0035] In formula (4), denotes the system state matrix at the th moment, denotes the system state matrix at the th moment, is a natural number greater than or equal to 1, denotes the state transition matrix at the th moment, denotes the control input matrix at the th moment, denotes the control quantity at the th moment, denotes process noise, with an expected value of 0;

[0036] S32. Determine the observation model for the target angle as:

[0037] (5)

[0038] In formula (5), denotes the observation value or measurement value of the system state matrix at the th moment, corresponding to the target angle measurement data after wild value elimination, that is, denotes the target angle measurement data after wild value elimination at the th moment, that is, the turntable azimuth angle measurement data after wild value elimination or the turntable elevation angle measurement data after wild value elimination at the th moment; denotes the state observation matrix, denotes the measurement noise at the th moment;

[0039] S33. Determine the prediction equation for the target angle as:

[0040] (6)

[0041] The prediction equation shown in formula (6) is obtained based on the state model in the step S31, wherein, denotes the prior state estimation value at the th moment, denotes the optimal estimation value at the th moment, that is, the target angle estimated true data corresponding to the th moment, denotes the prior estimation covariance at the th moment, denotes the posterior estimation covariance at the th moment, that is, the target angle estimation data after wild value elimination corresponding to the optimal estimation covariance of each moment, A T denotes a state transition matrix, denotes a system noise covariance matrix, and other parameters are the same as those in equation (4).

[0042] S34. Based on the iteration initial value, the target angle measurement data after the wild value is eliminated, and the prediction equation, the optimal estimation value of each moment is iteratively calculated to obtain the target angle estimation real data corresponding to each moment; wherein the iteration initial value includes a preset initial optimal estimation value and an initial posterior estimation covariance .

[0043] Further, the step S34 can include the following steps:

[0044] a1. determining a constant state transition matrix , a control input matrix , and a state observation matrix , determining a system noise covariance matrix , determining an observation noise covariance , and collectively serving as estimation parameters;

[0045] a2. For the target angle, determining the optimal estimation value and the posterior estimation covariance of the current iteration, which together constitute the input value of the current iteration; wherein the optimal estimation value and the posterior estimation covariance of the first iteration are respectively a preset initial optimal estimation value and an initial posterior estimation covariance ;

[0046] a3. Based on the estimation 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:

[0047] (7)

[0048] a4. Based on the calculated prior estimation covariance and the observation noise covariance , the Kalman gain is calculated using the Kalman gain calculation formula;

[0049] wherein the calculation formula of the Kalman gain is:

[0050] (8)

[0051] In equation (8), Indicates the first The Kalman gain at time step H T Represents the state observation matrix;

[0052] a5. Based on the calculated Kalman gain Prior state estimates and the Target angle measurement data after outlier removal at each time point The optimal estimate is calculated using the formula for calculating the optimal estimate. ;

[0053] The formula for calculating the optimal estimate is as follows:

[0054] (9)

[0055] In equation (9), the optimal estimate is Indicates the first The target angle at each moment is estimated using real data.

[0056] a6. Based on the calculated prior estimate covariance and Kalman gain Using the optimal estimation covariance calculation formula, the posterior estimation covariance is calculated. ;

[0057] The formula for calculating the optimal estimated covariance is as follows:

[0058] (10)

[0059] In equation (10), the posterior estimate of the covariance is... Indicates the first The optimal estimated covariance at each time step;

[0060] a7. The calculated optimal estimate and posterior estimate of covariance These together constitute the input values ​​for the next iteration. The iteration is repeated until the true target angle estimation data corresponding to each time step is calculated.

[0061] Furthermore, step S34 may also include the following steps:

[0062] b1. Determine the state transition matrix Control input matrix State observation matrix Determine the system noise covariance matrix based on the accuracy of the control system. The observation noise covariance is determined based on the observation model and sensor accuracy. , together serve as estimation parameters;

[0063] b2. For the target angle, determine the optimal estimate for the current iteration. and posterior estimate of covariance Together, they constitute the input value for the current iteration;

[0064] Among them, the optimal estimate corresponding to the first iteration and posterior estimate of covariance These are the preset initial optimal estimates. and initial posterior estimated covariance ;

[0065] b3. Based on the estimated parameters, the input value of the current iteration, and the updated prediction equation, calculate the prior state estimate of the target angle. and prior estimate of covariance ;

[0066] The updated prediction equation is as follows:

[0067] (11);

[0068] b4. Based on the calculated prior estimate covariance and the observed noise covariance The Kalman gain was calculated using the updated Kalman gain calculation formula. ;

[0069] The updated Kalman gain calculation formula is as follows:

[0070] (12);

[0071] b5. Based on the calculated Kalman gain Prior state estimates and the Target angle measurement data after outlier removal at each time point The optimal estimate is calculated using the updated optimal estimate calculation formula. , indicating the first Estimate the true data of the target angle at each moment;

[0072] The formula for calculating the updated optimal estimate is as follows:

[0073] (13);

[0074] b6. Based on the calculated prior estimate covariance and Kalman gain , the posterior estimation covariance is calculated by using the updated optimal estimation covariance calculation formula , representing the optimal estimation covariance at the th time point;

[0075] The updated optimal estimation covariance calculation formula is:

[0076] (14)

[0077] b7. The calculated optimal estimation value and the posterior estimation covariance together constitute the input value of the new iteration, and the iteration is repeated until the target angle estimation real data corresponding to each time point is calculated.

[0078] Further, the formula used in step S4 is:

[0079] (15)

[0080] In formula (15), represents the actual similarity calculated at the th time point, represents the azimuth angle theoretical data of the turntable at the th time point, represents the azimuth angle estimation real data of the turntable at the th time point, represents the pitch angle theoretical data of the turntable at the th time point, represents the pitch angle estimation real data of the turntable at the th time point, is a natural number.

[0081] Further, in step S5, according to all time points and the corresponding actual similarities respectively, the similarity function is fitted by using the least square method.

[0082] Further, in step S6, for each time point in the N consecutive new time points, steps S1-S4 are performed to obtain the actual similarity at the time point; steps S1-S5 are performed to obtain the estimation similarity at the time point; and then the absolute value of the difference between the actual similarity and the estimation similarity at the time point is calculated as the similarity difference value.

[0083] Further, in the step S7, the step of judging whether the difference of the similarity of the N continuous new time points meets the abnormal motion condition comprises: judging whether M of the difference of the similarity of the N continuous new time points are greater than the preset threshold, and if yes, determining that the abnormal motion condition is met; Wherein, M is a natural number greater than 0, and N / 2<M<N.

[0084] Compared with the prior art, the present application has the following advantages:

[0085] 1. According to the present application, a mathematical model is established according to the tracking of the azimuth and elevation angle of the turntable during the imaging process, the difference between the measured data and the theoretical data is evaluated to obtain the actual similarity, and a linear function fitting is performed to predict the estimated similarity of the new time, and the difference between the estimated similarity and the actual similarity of the new time is used to quickly find the abnormal motion of the space target, which can improve the timeliness of the research and judgment and provide support for data utilization.

[0086] 2. Compared with the traditional abnormal motion detection method based on the position error of the target in the camera field of view, the present application does not need to complete the astronomical positioning, and the method is simpler, and the target orbit abnormal motion can be quickly detected while the target is imaged, the timeliness is higher, and the resource integration and utilization is realized.

[0087] 3. The method of the present application fills the blank of using small-aperture telescopes on orbiting satellites to detect and identify the orbit abnormal motion of the target, and can quickly detect the orbit abnormal motion of the target based on imaging auxiliary data, effectively reduces the risk of space target collision caused by orbit maneuver, and can improve the safety of space target on-orbit operation. BRIEF DESCRIPTION OF DRAWINGS

[0088] Figure 1 A flowchart of a space target abnormal motion identification method based on imaging auxiliary data provided by the embodiment of the present application.

[0089] Figure 2 An iteration process diagram for obtaining target angle estimation real data of the embodiment of the present application.

[0090] Figure 3 Another iteration process diagram for obtaining target angle estimation real data of the embodiment of the present application.

[0091] Figure 4 A process understanding diagram of a space target abnormal motion identification method based on imaging auxiliary data provided by the embodiment of the present application. DETAILED DESCRIPTION

[0092] The technical solutions in the embodiments of the present application will be clearly and completely described with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, rather than all the embodiments of the present application. Based on the embodiments in the present application, all other embodiments obtained by a person of ordinary skill in the art without creative effort belong to the scope of protection of the present application.

[0093] Referring to Figure 1 , step 1: for each time, the azimuth angle theoretical data and the elevation angle theoretical data of the space target of a turntable are calculated according to a target orbit, and the azimuth angle measurement data and the elevation angle measurement data of the space target in a tracking imaging process of the turntable are obtained.

[0094] The space target in the embodiments of the present application can be various spacecrafts in space, which are not specifically limited here.

[0095] Specifically, step 1 includes: calculating the azimuth angle theoretical data and the elevation angle theoretical data of the space target at each time by using the orbit of the space target and the azimuth and elevation measurement values of the turntable, and obtaining the azimuth angle measurement data of the space target at each time output by the azimuth encoder of the turntable in real time and the elevation angle measurement data of the space target at each time output by the elevation encoder of the turntable in real time.

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

[0097] Through step 1, for each time, the azimuth angle theoretical data, the elevation angle theoretical data, the azimuth angle measurement data and the elevation angle measurement data of the time can be obtained, which belong to the imaging auxiliary data in the embodiments of the present application.

[0098] Step 2: the obtained azimuth angle measurement data and the obtained elevation angle measurement data are respectively subjected to outlier elimination by using an improved amplitude limiting filtering algorithm.

[0099] Since the obtained azimuth angle measurement data and the obtained elevation angle measurement data may contain outliers, that is, the outliers will have a great influence on the abnormal motion discrimination of the space target, the obtained azimuth angle measurement data and the obtained elevation angle measurement data are respectively subjected to outlier elimination by using an improved amplitude limiting filtering algorithm.

[0100] The embodiments of the present application first propose the following initial formula for outlier elimination on the basis of the amplitude limiting filtering algorithm:

[0101] (1)

[0102] In formula (1), indicates the data after outlier elimination, represents the sampling data corresponding to the th moment, represents the sampling data corresponding to the th moment, is a natural number greater than or equal to 1.

[0103] represents the threshold value corresponding to the sampling data, when the sampling data is the azimuth angle measurement data of the turntable, the data after the outlier elimination is the outlier-eliminated azimuth angle measurement data of the turntable, is the threshold value corresponding to the azimuth angle measurement data of the turntable; when the sampling data is the elevation angle measurement data of the turntable, the data after the outlier elimination is the outlier-eliminated elevation angle measurement data of the turntable, is the threshold value corresponding to the elevation angle measurement data of the turntable.

[0104] For the azimuth angle measurement data and the elevation angle measurement data of the turntable, the corresponding threshold values are not the same, for example, for the azimuth angle measurement data and the elevation angle measurement data of the turntable, the corresponding threshold values can be the maximum rotation angular velocity of the azimuth angle of the turntable and the maximum rotation angular velocity of the elevation angle of the turntable, respectively.

[0105] By setting , the outlier elimination can be performed on the azimuth angle measurement data and the elevation angle measurement data of the turntable corresponding to each current moment , respectively, to obtain the corresponding .

[0106] In order to set , assuming that the rotation angular acceleration of the azimuth angle of the turntable or the rotation angular acceleration of the elevation angle of the turntable is , then:

[0107] (2)

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

[0109] Then, the improved amplitude limiting filtering algorithm is used for the outlier elimination of the obtained azimuth angle measurement data and the obtained elevation angle measurement data of the turntable, and the formula used can be:

[0110] (3)

[0111] In formula (3), represents the sampling data corresponding to the The sampling data corresponding to each time point, and the other parameters are the same as in equation (1).

[0112] Traditional methods do not preprocess the turntable azimuth and pitch angle measurement data, which are used as raw data. Outliers in these raw data may lead to misjudgments of space target anomalies. At the same time, the error between the measured value and the true value may prevent the judgment from being completed. This invention proposes an improved amplitude limiting filtering algorithm to preprocess the raw data by removing outliers. This can improve the efficiency and accuracy of data processing, lay a solid foundation for subsequent data processing, effectively avoid misjudgments caused by outliers, and improve the accuracy of space target anomaly detection.

[0113] Step 3: Filter and estimate the turntable azimuth and pitch measurement data after outlier removal to obtain the actual estimated data of the turntable azimuth and pitch at each time point.

[0114] In one optional implementation, step 3 may include steps 31 to 34.

[0115] In this embodiment of the invention, steps 31 to 34 are executed for the turntable azimuth angle and the turntable pitch angle, respectively, to obtain the actual estimated data of the turntable azimuth angle and the actual estimated data of the turntable pitch angle at each time.

[0116] In steps 31 to 34, the target angle is described as either the turntable azimuth angle or the turntable pitch angle.

[0117] Step 31: Determine the state model for the target angle, i.e., the turntable azimuth or turntable pitch angle, as follows:

[0118] (4)

[0119] In equation (4), Indicates the first The system state matrix at each time step Indicates the first The system state matrix at each time step A natural number greater than or equal to 1. Indicates the first The state transition matrix at each time step Indicates the first The control input matrix at each time step Indicates the first Control quantity at any given moment This represents process noise, with an expected value of 0.

[0120] Step 32: Determine the observation model for the target angle as follows:

[0121] (5)

[0122] In equation (5), Indicates the first The observed or measured values ​​of the system state matrix at each moment correspond to the target angle measurement data after outlier removal. In other words... Indicates the first The target angle measurement data after outlier removal at time step i, i.e., the i-th time step i. The turntable azimuth measurement data or the turntable pitch measurement data after outlier removal at any given moment; Represents the state observation matrix, Indicates the first Measurement noise at each moment.

[0123] Step 33: Determine the prediction equation for the target angle as follows:

[0124] (6)

[0125] The prediction equation shown in equation (6) is obtained based on the state model in step S31, wherein, Indicates the first The prior state estimate at time 1. Indicates the first The optimal estimate at time i, i.e., the i-th time... The actual data for the target angle estimation at each time point. Indicates the first The prior estimate of covariance at time 1. Indicates the first The posterior estimated covariance at time i, i.e., the i-th time... The optimal estimated covariance at time A T Represents the state transition matrix. Let represent the system noise covariance matrix, and the other parameters are the same as in equation (4).

[0126] Step 34: Based on the initial iteration value, the target angle measurement data after outlier removal, and the prediction equation, iteratively calculate the optimal estimate value at each time step to obtain the true target angle estimate data corresponding to each time step; wherein, the initial iteration value includes a preset initial optimal estimate value. and initial posterior estimated covariance .

[0127] For details, see Figure 2 Step S34 may include the following steps:

[0128] Step a1: Determine the constant state transition matrix , control input matrix , and state observation matrix , determine system noise covariance matrix , determine observation noise covariance , together as an estimated parameter.

[0129] Step a2: for the target angle, determine the optimal estimated value corresponding to the current iteration and the posteriori estimated covariance , together constitute the input value of the current iteration; wherein the optimal estimated value corresponding to the first iteration and the posteriori estimated covariance are respectively preset initial optimal estimated value and initial posteriori estimated covariance .

[0130] Step a3: based on the estimated parameter, the input value of the current iteration and the prediction equation, calculate the prior state estimated value and the prior estimated covariance of the target angle, to obtain the following results:

[0131] (7).

[0132] Step a4: based on the calculated prior estimated covariance and the observation noise covariance , calculate the Kalman gain using the Kalman gain calculation formula;

[0133] Wherein, the calculation formula of the Kalman gain is:

[0134] (8)

[0135] In formula (8), represents the Kalman gain at the th moment, H T represents the state observation matrix.

[0136] Step a5: based on the calculated Kalman gain , the prior state estimated value and the target angle measurement data after wild value elimination at the th moment , calculate the optimal estimated value using the optimal estimated value calculation formula;

[0137] Wherein, the optimal estimated value calculation formula is:

[0138] (9)

[0139] in formula (9), the optimal estimation value represents the target angle estimation real data at the th moment.

[0140] Step a6: based on the calculated prior estimation covariance and the Kalman gain , the posterior estimation covariance is calculated by using the optimal estimation covariance calculation formula.

[0141] The optimal estimation covariance calculation formula is as follows:

[0142] (10)

[0143] In formula (10), the posterior estimation covariance represents the optimal estimation covariance at the th moment.

[0144] It can be understood that, for the th moment, the optimal estimation value can be obtained through the above steps a1~a6, as the target angle estimation real data at the th moment, and the posterior estimation covariance is obtained, which is equivalent to completing the real data estimation at the th moment.

[0145] Step a7: the calculated optimal estimation value and the posterior estimation covariance constitute the input value of a new iteration, and the iteration is repeated until the target angle estimation real data corresponding to each moment is calculated.

[0146] The calculated optimal estimation value and the posterior estimation covariance constitute the input value of a new iteration, and the steps a2~a6 are executed, the target angle estimation real data at the next moment can be obtained, and the posterior estimation covariance is obtained, which is equivalent to completing the real data estimation at the th moment.

[0147] Then, through continuous iteration, the azimuth estimation real data and the pitch estimation real data of the turntable corresponding to each moment in step 1 can be obtained.

[0148] The embodiment of the application can select a suitable initial optimal estimation value and the initial posteriori estimation covariance As the iterative initial value, the iterative initial value and the estimation parameter are substituted into the filtering equation constituted by formula (7) to formula (10) to carry out iteration, so that the measurement value is more close to the true value. It can be understood that the above filtering equation is realized based on Kalman filtering.

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

[0150] Wherein, the ideal physical model of the angle of the azimuth angle of the turntable and the pitch angle of the turntable is as follows:

[0151] (11)

[0152] In formula (11), indicates the azimuth angle of the turntable or the pitch angle of the turntable at the i th moment, indicates the azimuth angle of the turntable or the pitch angle of the turntable at the i th moment, indicates the angular velocity of the azimuth angle of the turntable or the pitch angle of the turntable at the i th moment, indicates the angular acceleration, indicates the sampling period.

[0153] Convert formula (11) into a matrix form, and then there is:

[0154] (12)

[0155] Let , , then there is:

[0156] (13)

[0157] In the actual system, the satellite remote sensing data only contains angle information, and does not contain angular velocity and angular acceleration information, so let the angular velocity , and the angular acceleration , so formula (12) can be simplified as:

[0158] (14)

[0159] In formula (14), indicates a selected value, and the value range is .

[0160] Therefore, it can be concluded that in the actual system: the state transition matrix , the control input matrix , and the system noise covariance matrix ​The specific values are determined according to the control system accuracy. The acquired parameters are still the azimuth angle of the turntable and the pitch angle of the turntable, and therefore, the observation model corresponding to formula (5) can be changed to:

[0161] (15)

[0162] wherein, represents a selected value, and the value range is

[0163] Therefore, in the actual system, the state observation matrix , and the observation noise covariance The specific values are determined according to the observation model and the sensor accuracy. Therefore, compared with steps a1 to a7, the filtering equation can be simplified for the actual system.

[0164] Specifically, referring to Figure 3 , step 34 can include the following steps b1 to b7:

[0165] Step b1: determining a state transition matrix , a control input matrix , a state observation matrix , and a system noise covariance matrix according to the control system accuracy, and an observation noise covariance according to the observation model and the sensor accuracy, which are collectively used as estimation parameters.

[0166] Step b2: determining, for the target angle, an optimal estimation value and a posteriori estimation covariance of the current iteration, which collectively constitute the input value of the current iteration.

[0167] wherein, the optimal estimation value and the posteriori estimation covariance of the first iteration are respectively a preset initial optimal estimation value and an initial posteriori estimation covariance .

[0168] Step b3: based on the estimation parameters, the input value of the current iteration, and an updated prediction equation, calculating a priori state estimation value and a priori estimation covariance of the target angle.

[0169] wherein, the updated prediction equation is:

[0170] (16).

[0171] Step b4: Based on the calculated prior estimate of covariance and the observed noise covariance The Kalman gain was calculated using the updated Kalman gain calculation formula. ;

[0172] The updated Kalman gain calculation formula is as follows:

[0173] (17).

[0174] Step b5: Based on the calculated Kalman gain Prior state estimates and the Target angle measurement data after outlier removal at each time point The optimal estimate is calculated using the updated optimal estimate calculation formula. , indicating the first Estimate the true data of the target angle at each moment;

[0175] The formula for calculating the updated optimal estimate is as follows:

[0176] (18).

[0177] Step b6: Based on the calculated prior estimate of covariance and Kalman gain Using the updated optimal estimate covariance calculation formula, the posterior estimate covariance is calculated. , indicating the first The optimal estimated covariance at each time step;

[0178] The formula for calculating the updated optimal estimated covariance is as follows:

[0179] (19).

[0180] Understandably, regarding the first At any given time, the optimal estimate can be obtained through steps b1 to b6 described above. As the first Real data of target angle estimation at each moment And obtain the posterior estimate of covariance. Indicates the first The optimal covariance estimate at time step t was completed, thus achieving the t-th time step. Estimate the actual data at each moment.

[0181] Step b7: Calculate the optimal estimate. and the posterior estimation covariance The input values of the new iteration are composed of the estimated real data of the target angle and the posterior estimation covariance

[0182] The optimal estimation value calculated and the posterior estimation covariance The estimated real data of the target angle at the next moment and the posterior estimation covariance are obtained by executing steps b2 to b6 with the input values of the new iteration , which is equivalent to completing the real data estimation at the moment.

[0183] Then, through continuous iteration, the estimated real data of the azimuth angle of the turntable and the estimated real data of the elevation angle of the turntable at each moment corresponding to step 1 can be obtained.

[0184] In the embodiment of the application, the state transition matrix is determined, the control input matrix is determined, the state observation matrix is determined, the system noise covariance matrix is determined according to the control system accuracy, the observation noise covariance is determined according to the observation model and the sensor accuracy, the estimated parameters conforming to the actual system can be obtained, the estimated real data of the azimuth angle of the turntable and the estimated real data of the elevation angle of the turntable can be made closer to the real data, and thus the accuracy of subsequent space target anomaly discrimination can be improved.

[0185] In the embodiment of the application, the real values of the azimuth angle and the elevation angle of the turntable can be effectively estimated based on the Kalman filtering algorithm by inputting the measurement values, i.e., the measurement data of the azimuth angle of the turntable and the measurement data of the elevation angle of the turntable, the measurement values can be corrected, and the relative error between the measurement values and the real values can be reduced.

[0186] Step 4: The actual similarity between the estimated real data and the theoretical data at the corresponding moment is calculated by using the obtained theoretical data of the azimuth angle of the turntable, the estimated real data of the azimuth angle of the turntable, the theoretical data of the elevation angle of the turntable and the estimated real data of the elevation angle of the turntable at each moment.

[0187] The formula used in step 4 is:

[0188] (20)

[0189] In formula (20), represents the actual similarity calculated at the moment, represents the theoretical data of the azimuth angle of the turntable at the moment, represents the azimuth angle estimation real data of the turntable at the th moment, represents the pitch angle theoretical data of the turntable at the th moment, represents the pitch angle estimation real data of the turntable at the th moment, is a natural number.

[0190] It can be understood that formula (20) is realized based on the Euclidean distance calculation formula. For each moment , the corresponding actual similarity can be calculated, and thus a plurality of sets of data can be obtained, wherein represents the moment .

[0191] Step 5: fitting a similarity function according to the calculated actual similarities of all moments;

[0192] In an optional embodiment, step 5 can include:

[0193] fitting a similarity function by using the least square method according to all moments and the respective corresponding actual similarities, for predicting the estimation similarity of a new moment.

[0194] Specifically, the plurality of sets of data obtained in step 4 can be fitted into a linear equation representing the similarity function by using the least square method, which can be specifically:

[0195] (21)

[0196] In formula (21), represents the estimation similarity, that is, the similarity between the estimation real data and the theoretical data of the th moment predicted by inputting a moment into the similarity function, and are fitting coefficients, which can be determined by the least square method.

[0197] The embodiment of the application fits a similarity function by using the actual similarities of a plurality of historical moments, can predict the corresponding estimation similarity for each new moment, and can subsequently determine the abnormal situation of the space target by using the difference between the estimation similarity of a plurality of new moments and the actual similarity calculated by the real.

[0198] Step 6: obtaining the actual similarity of each of the N continuous new time points, obtaining the estimated similarity of the time point by using the similarity function, and calculating the absolute value of the difference between the actual similarity and the estimated similarity at the time point as the similarity difference value;

[0199] N is a natural number greater than 0, which can be selected as required, such as 8, 10, etc.

[0200] It should be noted that the N continuous new time points in step 6 refer to time points that have not been calculated for actual similarity, have not been used to fit the similarity function, and have not been used to determine spatial target anomalies. For the first time, step 6 is executed, the N continuous new time points refer to the multiple groups of data used to fit the similarity function N continuous new time points after that.

[0201] For each of the N continuous new time points, the actual similarity of the time point can be obtained according to steps 1-4, which will not be repeated here.

[0202] For each of the N continuous new time points, the time point can be substituted into the similarity function shown in equation (21) to obtain the estimated similarity of the time point.

[0203] Then, the absolute value of the difference between the actual similarity and the estimated similarity at the time point is calculated as the similarity difference value.

[0204] Step 7: determining whether the similarity difference values of the N continuous new time points meet the anomaly condition, if so, determining that the spatial target has an anomaly; if not, returning to step 6.

[0205] In an optional embodiment, the determination of whether the similarity difference values of the N continuous new time points meet the anomaly condition includes: determining whether M of the similarity difference values of the N continuous new time points are all greater than a preset threshold, if so, determining that the anomaly condition is met; wherein M is a natural number greater than 0, and N / 2

[0206] It can be understood that the anomaly condition is met when more than half of the similarity difference values of the N continuous new time points are greater than the preset threshold, which can avoid misjudgment.

[0207] In the specific implementation process, each of the N continuous new time points can be sequentially determined, if the similarity difference value of the current time point is greater than the preset threshold, the counter is incremented by one; when the counter value reaches M, it is determined that the anomaly condition is met, and it is determined that the spatial target has an anomaly; otherwise, return to step 6 to continue monitoring and anomaly determination for each of the N continuous new time points that appear again.

[0208] As a preferred implementation, N can be 10, and M can be 6, that is, 10 is judged for the similarity difference value calculation result, and 6 is judged, that is, if 6 of the 10 consecutive groups of similarity difference values are greater than the preset threshold value, it is determined that the space target has moved, otherwise, step 6 is returned to continue monitoring and discriminating.

[0209] The implementation process of the space target movement discrimination method based on imaging auxiliary data provided by the embodiment of the application is described in detail below with reference to Figure 4 , and the specific content is described above, which is not repeated here.

[0210] The embodiment of the application fits a similarity function through a large amount of data of the space target that has not changed the orbit, and uses the similarity function to evaluate whether the space target has changed the orbit through state estimation, and specifically, the corresponding similarity difference value is obtained by comparing the actual similarity and the estimated similarity at each time, and compared with the preset threshold value, when the similarity difference value is greater than the preset threshold value at a certain number of times, it is determined that the space target has moved. This method does not need to complete astronomical positioning, and the method is simpler, and the target orbit movement can be quickly detected while the target is imaged, and the timeliness is higher.

[0211] In the space target movement discrimination method based on imaging auxiliary data provided by the embodiment of the application, for the azimuth and elevation angles of the space target, first, the theoretical data and the measured data at each time are obtained; then the measured data is subjected to wild value elimination using an improved amplitude limiting filtering algorithm to improve the data processing efficiency and accuracy, and to avoid misjudgment caused by wild values; next, the measured data after wild value elimination is subjected to filtering estimation to obtain the estimated real data corresponding to each time to correct the measured data and reduce the relative error; then the actual similarity between the estimated real data and the theoretical data at each time is calculated; and a similarity function is fitted according to the calculated actual similarity at all times; then for each time in N consecutive new times, the actual similarity at the time is obtained, the estimated similarity at the time is obtained using the similarity function, the absolute value of the difference between the actual similarity and the estimated similarity at the time is calculated as a similarity difference value, and whether the similarity difference values of the N consecutive new times meet the movement condition is judged to discriminate whether the space target has moved.

[0212] The embodiment of the present application is according to tracking the azimuth and elevation angle change of the turntable in the imaging process, establishing a mathematical model, evaluating the gap between the measured data and the theoretical data to obtain the actual similarity, and performing linear function fitting to predict the estimated similarity at a new time, and using the difference between the estimated similarity at the new time and the actual similarity to quickly find the abnormal movement of the space target. The method can improve the timeliness of research and judgment, and provide support for data application. Moreover, compared with the traditional abnormal movement detection method based on the position error of the target in the camera field of view, the present application does not need to complete astronomical positioning, and the method is simpler. The target track abnormal movement can be quickly detected while the target is being imaged, the timeliness is higher, and resource integration and application are realized.

[0213] The method of the present application fills the blank of detecting and discriminating the track abnormal movement of the target by using the small-aperture telescope of the on-orbit satellite, can quickly detect the track abnormal movement of the target based on the imaging auxiliary data, effectively reduces the space target collision risk caused by the track maneuver, and can improve the safety of the on-orbit operation of the space target.

[0214] Although the present specification is described in terms of embodiments, not every embodiment contains only one independent technical solution. The description of the specification is only for the sake of clarity, and those skilled in the art should consider the specification as a whole. The technical solutions in each embodiment can also be appropriately combined to form other embodiments that those skilled in the art can understand.

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

Claims

1. A method for identifying a spatial target anomaly based on imaging auxiliary data, characterized in that: Specifically comprising the following steps: S1. For each time, calculate the azimuth angle theoretical data and the elevation angle theoretical data of the space target according to the target orbit, and obtain the azimuth angle measurement data and the elevation angle measurement data of the space target in the tracking imaging process; S2. The obtained azimuth angle measurement data and the elevation angle measurement data are respectively subjected to wild value elimination by using the improved amplitude limiting filtering algorithm; S3. The azimuth angle measurement data and the elevation angle measurement data after wild value elimination are respectively subjected to filtering estimation, so as to obtain the estimated real data of the azimuth angle and the estimated real data of the elevation angle corresponding to each time; S4. The actual similarity between the estimated real data and the theoretical data at the corresponding time is calculated by using the obtained azimuth angle theoretical data, the estimated real data of the azimuth angle, the elevation angle theoretical data and the estimated real data of the elevation angle at each time; S5. The similarity function is fitted according to the actual similarity at all times, and the estimated similarity at the time is calculated; S6. For each time in N continuous new times, the actual similarity at the time is obtained, the estimated similarity at the time is obtained by using the similarity function, and the absolute value of the difference between the actual similarity and the estimated similarity at the time is calculated as the similarity difference value, wherein N is a natural number greater than 0; S7. It is judged whether the similarity difference values of the N continuous new times meet the abnormal motion condition, if yes, it is determined that the space target has an abnormal motion, and if not, it is returned to step S6.

2. The method of claim 1, wherein the method comprises: The step S1 specifically comprises: calculating the azimuth angle theoretical data and the elevation angle theoretical data of the space target at each time by using the orbit of the space target and the azimuth and elevation measurement values of the turntable, and obtaining the azimuth angle measurement data of the space target at each time output by the azimuth encoder of the turntable in real time, and the azimuth angle measurement data of the space target at each time output by the elevation encoder of the turntable in real time.

3. The method of claim 1, wherein the method further comprises: determining a plurality of candidate target objects in the image based on the auxiliary data; and determining a target object in the image based on the plurality of candidate target objects. The step S2 specifically comprises the following steps: S21. The following initial formula for wild value elimination is proposed on the basis of the amplitude limiting filtering algorithm: （1） In formula (1), represents the data after wild value elimination, represents the sampling data corresponding to the time point, represents the sampling data corresponding to the time point, is a natural number greater than or equal to 1. represents a threshold value corresponding to the sampling data, when the sampling data is the azimuth angle measurement data of the turntable, the data after the wild value elimination is the wild value eliminated azimuth angle measurement data of the turntable, represents a threshold value corresponding to the sampling data; when the sampling data is the pitch angle measurement data of the turntable, the data after the wild value elimination is the wild value eliminated pitch angle measurement data of the turntable, represents a threshold value corresponding to the sampling data; when the sampling data is the pitch angle measurement data of the turntable, the data after the wild value elimination is the wild value eliminated pitch angle measurement data of the turntable, S22. Set the value of the variable θ in equation (1) to : Assuming that the rotational angular acceleration of the azimuth angle of the turret or the elevation angle of the turret is , then: （2） In formula (2), denotes a sampling interval, y n-1 denotes the sampling data corresponding to the (n-1)th moment, y n-2 denotes the sampling data corresponding to the (n-2)th moment; S23. Combined with formula (1) and (2), the azimuth angle measurement data and the elevation angle measurement data obtained in step S1 are respectively subjected to wild value elimination by using the improved amplitude limiting filtering algorithm, and the formula is as follows: （3） In formula (3), represents the sampling data corresponding to the first moment, and the remaining parameters are the same as those in formula (1).

4. The method of claim 1, wherein the method further comprises: The step S3 specifically comprises the following steps: S31. It is determined that the state model for the target angle, i.e. the azimuth angle or the elevation angle of the turntable, is: （4） In formula (4), denotes the system state matrix at the time point, denotes the system state matrix at the time point, is a natural number greater than or equal to 1, denotes the state transition matrix at the time point, denotes the control input matrix at the time point, denotes the control quantity at the time point, denotes the process noise, and the expected value is 0; S32. It is determined that the observation model for the target angle is: （5） In formula (5), represents the observation or measurement value of the system state matrix at the time, corresponding to the target angle measurement data after outlier rejection, that is, represents the target angle measurement data after outlier rejection at the time, that is, the target angle measurement data after outlier rejection at the time or the target elevation angle measurement data after outlier rejection at the represents the state observation matrix, represents the measurement noise at the time. S33. It is determined that the prediction equation for the target angle is: （6） The prediction equation shown in formula (6) is based on the state model in step S31, wherein, represents the prior state estimation value at the th moment, represents the optimal estimation value at the th moment, that is, the target angle estimation real data corresponding to the th moment, represents the prior estimation covariance at the th moment, represents the posterior estimation covariance at the th moment, that is, the optimal estimation covariance at the th moment, A T represents the state transition matrix, represents the system noise covariance matrix, and the remaining parameters are the same as formula (4); S34. Based on the iteration initial value, the target angle measurement data after the wild value is eliminated, and the prediction equation, the optimal estimation value at each time is iteratively calculated to obtain the target angle estimation real data corresponding to each time; wherein the iteration initial value includes a preset initial optimal estimation value and an initial posterior estimation covariance .

5. The method of claim 4, wherein the method further comprises: The step S34 can comprise the following steps: a1. determining a constant state transition matrix , a control input matrix , and a state observation matrix , determining a system noise covariance matrix , determining an observation noise covariance , jointly as estimation parameters; a2. determining an optimal estimation value corresponding to the current iteration for the target angle and a posterior estimation covariance together constitute an input value of the current iteration; wherein the optimal estimation value and the posterior estimation covariance corresponding to the first iteration are respectively a preset initial optimal estimation value and an initial posterior estimation covariance ; a3. Based on the estimated parameter, 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, and the following results are obtained: （7） a4. based on the computed a priori estimate covariance and the observation noise covariance , the Kalman gain is computed using the Kalman gain computation formula ; Wherein, the calculation formula of the Kalman gain is: （8） In formula (8), denotes the Kalman gain at the time point, H T denotes the state observation matrix; a5. based on the calculated Kalman gain , the prior state estimate , and the outlier-removed target angle measurement data at the first time , the optimal estimate is calculated using the optimal estimate calculation formula ; Wherein, the calculation formula of the optimal estimation value is: （9） In formula (9), the optimal estimation value represents the target angle estimation real data at the first time point; a6. Based on the computed prior estimate covariance and the Kalman gain , the posterior estimate covariance is computed using the optimal estimate covariance formula ; Wherein, the calculation formula of the optimal estimation covariance is: （10） In Equation (10), the posteriori estimation covariance represents the optimal estimation covariance at the th time instant; a7. The computed optimal estimate and the posteriori estimate covariance together constitute the input values for a new iteration, which is repeated until the target angle estimate for each time instant is computed.

6. The method of claim 4, wherein the method further comprises: The step S34 can also comprise the following steps: b1. determining a state transition matrix , a control input matrix , a state observation matrix , determining a system noise covariance matrix based on control system accuracy , determining an observation noise covariance based on the observation model and sensor accuracy , jointly as estimation parameters; b2. determining an optimal estimation value corresponding to the current iteration for the target angle and a posteriori estimation covariance together constitute an input value of the current iteration; Wherein, the optimal estimation value corresponding to the first iteration and the posterior estimation covariance are respectively a preset initial optimal estimation value and an initial posterior estimation covariance ; b3. based on the estimated parameters, the input values of the current iteration, and the updated prediction equation, compute a priori state estimate values of the target angle and a priori estimate covariance ; Wherein, the updated prediction equation is: （11）； b4. Based on the computed a priori estimate covariance and the observation noise covariance , the Kalman gain is computed using the updated Kalman gain computation formula ; Wherein, the calculation formula of the updated Kalman gain is: （12）； b5. based on the calculated Kalman gain , the prior state estimate , and the outlier-removed target angle measurement data at the first time , the optimal estimate is calculated using the updated optimal estimate calculation formula , which represents the true data of the target angle estimate at the first time The updated optimal estimation value calculation formula is: （13）； b6. Based on the computed a priori estimation covariance and the Kalman gain , the a posteriori estimation covariance is computed using the updated optimal estimation covariance computation formula , which represents the optimal estimation covariance at the th time instant; The updated optimal estimation covariance calculation formula is: （14）； b7. The computed optimal estimate and the posteriori estimate covariance together constitute the input values for a new iteration, which is repeated until the target angle estimate for each time instant is computed.

7. The method of claim 1, wherein the method further comprises: determining a plurality of candidate objects in the image based on the auxiliary data; and determining a plurality of candidate objects in the image based on the auxiliary data. The formula used in step S4 is: （15） In equation (15), Indicates the first The actual similarity calculated at each time point. Indicates the first Theoretical data of the turntable azimuth angle at each moment Indicates the first The estimated true data of the turntable azimuth angle at each moment. Indicates the first Theoretical data of the turntable pitch angle at a given moment Indicates the first The estimated true data of the turntable pitch angle at each moment. It is a natural number.

8. The method of claim 1, wherein the method further comprises: determining a plurality of candidate objects in the image; and determining a plurality of candidate object trajectories based on the plurality of candidate objects. In step S5, the similarity function is fitted by using the least square method according to all time points and the corresponding actual similarities.

9. The method of claim 1, wherein the method further comprises: determining a plurality of candidate objects in the image based on the auxiliary data; and determining a plurality of candidate objects in the image based on the auxiliary data. In step S6, for each time point in the N continuous new time points, steps S1-S4 are performed to obtain the actual similarity of the time point, and steps S1-S5 are performed to obtain the estimated similarity of the time point; then the absolute value of the difference between the actual similarity and the estimated similarity at the time point is calculated as the similarity difference value.

10. The method of claim 1, wherein the method further comprises: determining a plurality of candidate objects in the image; and determining a plurality of candidate object trajectories based on the plurality of candidate objects. In step S7, it is judged whether the similarity difference values of the N continuous new time points satisfy the abnormal condition, including: judging whether M similarity difference values in the N continuous new time points are greater than a preset threshold, if yes, it is determined that the abnormal condition is satisfied; wherein M is a natural number greater than 0, and N / 2

Images