An Improved Laplace Satellite Tracking Method Based on the Kalman Filter

Abstract

When photoelectric measuring equipment is used to track satellites, the extraction of the short-term or long-term target often fails because the target is weak, clouds block the target, and/or the sun’s angle is too small, resulting in the loss of the tracking target. In this study, an improved Laplacian satellite tracking method based on the Kalman filter is proposed. Firstly, the improved Laplacian algorithm was used for the initial fitting of the equation of motion of a small amount of measurement data. Judgment of the validity and Kalman filtering was carried out on the current frame’s measurement data to calculate the optimal estimate of the current frame’s orbit data, and the accurate equation of motion was iteratively fitted to obtain high-precision data for predicting the satellite’s orbit frame by frame. Numerical tracking of the equipment was carried out. This method was experimentally validated on an actual optical measurement device. The test results showed that this method can make up for the frequent loss of short-term targets. Under the condition that the maximum deviation is less than 3″, the length of extrapolated data can be up to 30 s and the length of the measurement data was less than 30 s. This method may improve the stability of tracking equipment as well as the accuracy and integrity of the measurement data.

1. Introduction

With the progress of science and technology, satellites have become an important means of obtaining all kinds of information; they are the main source of collecting important information and are favored by all countries. As an effective means of directly obtaining high-resolution image information, imaging satellites can carry out all-weather, all-day, and high-resolution real-time reconnaissance on various strategic and tactical targets and are the main source of important intelligence [1,2,3]. At the same time, the tracking, measurement, maintenance, and control of satellites have become the focus of new research [4,5,6]. Photoelectric measurement equipment has become an important tool for tracking and measuring satellites because of its capability for high-precision tracking and measurement [7,8,9].

Horizontal photoelectric measuring equipment is usually the first choice for tracking and measuring artificial satellites because it can avoid the zenith’s blind area. Perturbations Satellite Orbit Model 4 (SGP4) [10,11] is used to calculate a satellite’s orbit according to the recent bivariate elements of the orbit [12,13,14] obtained by the observatory. The calculated orbit data are used to guide the servo control device and drive the observation device to search for and capture the satellite. After the satellite enters the field of view of the observation camera, the image processing device extracts the miss distance of the target in the field of view, which is used as the tracking error input of the servo control device, and switches the automatic control closed-loop tracking to start accurate tracking of the satellite. The accuracy of the miss distance for that target extracted by the image processing equipment determines the tracking ability of the automatic closed-loop servo control. Several consecutive frames of an image or a long-term target’s miss distance with an extraction error will lead to the loss of the tracking target. In the process of tracking small, dark, low-orbit satellites, the target is often dark, and the signal-to-noise ratio of successive frames of images is too low, which leads to the situation that the short-term image processing equipment cannot accurately extract the target’s miss distance within tens of milliseconds. It may also occur that the target is blocked by clouds or the sun’s angle is too small, meaning that the long-term image processing equipment cannot accurately extract the target’s miss distance within tens of seconds, eventually leading to the loss of the tracking target, which seriously affects the stability and integrity of the satellite tracking and measurement system. In this case, if the measurement data of the previous section are used for extrapolation of the track to predict the subsequent data of the track and realize high-precision data for guiding the servo control equipment, this can avoid the situation in which the short-term target’s miss distance cannot be accurately extracted because of the unstable extraction of the weak target. In the case where the long-term target’s miss distance cannot be accurately extracted due to cloud occlusion or the sun’s angle being too small, long-term accurate extrapolation data can be provided. At this time, the target has not been lost and is still near the center of the field of view. When the target appears in the field of view again and the image processing equipment can extract the target’s miss distance, accurate tracking of the target can be continued [10,11]. However, these methods in the literature have either used multi-station data, which are unable to meet the needs of a single station, or need to use long-term tracking data, and the accuracy is low. This study focuses on single-station equipment for accurate predictions over a short arc.

For predicting the orbit by using single-station photoelectric equipment, the improved Laplacian algorithm has mostly been used in the literature [15,16]. The main advantage of this method is that extrapolation of the orbit can be carried out by using the precise tracking data obtained by single-station photoelectric measuring equipment, and information on the distance between the satellite and the measuring station is not required. This overcomes the difficulty of photoelectric equipment in obtaining information on the satellite’s distance.

When the improved Laplacian algorithm is used for predicting orbits, a longer period of tracking and measurement data is required as a prerequisite to extract long-term and high-precision tracking data. The main reason is that in the process of tracking and measurement, the deviation in the extraction of faint targets and other aspects creates noise with a large amplitude, that is, the accuracy of the measurement data is low. Noise will affect the accuracy of the fit of the equation of motion and then affect the accuracy of the extrapolated tracking data. Data can be extrapolated with high accuracy only by fitting the equation of motion with the measurement data for a long time.

However, in actual work, due to the short transit time of low-orbit satellites, the observable time is generally only a few minutes, and it is difficult to meet the requirement for long-term tracking and measurement data in many cases. Therefore, this study proposed an improved Laplacian satellite tracking method based on the Kalman filter. On the one hand, the commonly used method of image processing via the automatic closed-loop tracking mode of the servo controller was cancelled. Instead, precise predicted data of the trajectory were used to guide the servo controller equipment throughout the process. On the other hand, the Kalman filter [17,18,19] was used to obtain the optimal estimation of the measurement data of tracking to improve their accuracy. The method was combined with the improved Laplacian algorithm to accurately fit the equation of motion and then obtain high-precision prediction data for the satellite’s orbit frame by frame. This can reduce the length of the measurement data used and improve the accuracy and time length of the prediction data.

We implemented this method on an actual optical measurement device and conducted experimental verification. The experimental results showed that this method can solve the problem of frequent target loss for a short time caused by the unstable extraction of the miss distance of weak targets and can significantly improve the stability of tracking. At the same time, in the case of long-term target loss, as long as a certain length of actual tracking data is obtained in advance and the optimal estimation of tracking data is obtained through the validity judgment of measurement data and the Kalman filter, high-precision data on the long-term orbit can be predicted. When the target appears in the field of view again, the method can still ensure that the tracking measurement continues near the center of the field of view and maximize the tracking and measurement work to provide more complete and reliable measurement data.

2. Materials and Methods

2.1. An Improved Laplace Algorithm

The improved Laplacian algorithm has been applied previously for predicting the orbit of satellites. The main advantage of this method is that it can use only the short-arc tracking data obtained by the photoelectric equipment of a single station to fit the equation of motion and carry out extrapolation of the orbit. Moreover, it does not need to know the distance between the satellite and the measuring station, which overcomes the problem that the photoelectric equipment cannot obtain information on the distance of the satellite. The main steps of this method are as follows.

Under the action of various forces of perturbation, the orbit of a satellite will have long-term, long-period, and short-period changes and even irregular changes. Different types of orbits are subject to different forces of perturbation. Low- and medium-orbit satellites are mainly affected by Earth’s non-spherical gravity and atmospheric drag, while for short-arc tracking measurements, the influence of long-period atmospheric drag can be ignored. Therefore, only the model of the satellite’s dynamics accounting for the Earth’s non-spherical gravity needs to be considered.

The potential function of the earth’s gravitational field is shown in Equation (1):

$$\mathrm{V}={\displaystyle \frac{\mathrm{GM}}{\mathrm{R}}}\left[1-{\displaystyle \sum }_{\mathrm{l}=2}^{\infty }{\mathrm{J}}_{\mathrm{l}}{\left({\displaystyle \frac{{\mathrm{a}}_{\mathrm{e}}}{\mathrm{R}}}\right)}^{\mathrm{l}}{\mathrm{P}}_{\mathrm{l}}\left(\sin \mathsf{\phi }\right)+{\displaystyle \sum }_{\mathrm{l}=2}^{\infty }{\displaystyle \sum }_{\mathrm{m}=1}^{\mathrm{l}}{\left({\displaystyle \frac{{\mathrm{a}}_{\mathrm{e}}}{\mathrm{R}}}\right)}^{\mathrm{l}}{\mathrm{P}}_{\mathrm{lm}}\left(\sin \mathsf{\phi }\right)\left({\mathrm{C}}_{\mathrm{lm}}\cos \mathrm{m}\mathsf{\lambda }+{\mathrm{S}}_{\mathrm{lm}}\sin \mathrm{m}\mathsf{\lambda }\right)\right]$$

(1)

${\mathrm{P}}_{\mathrm{lm}}$ is the associated Legendre function. ${\mathrm{C}}_{\mathrm{lm}}$ and ${\mathrm{S}}_{\mathrm{lm}}$ are harmonic coefficients of non-normalized gravitational potential expansions. When m = 0, it is transformed into a zonal harmonic term, and when m ≠ 0, it is called a field harmonic term. For the orbit prediction of a low-orbit satellite, only the harmonic term corresponding to m = 0 can meet the accuracy requirements.

The gravitational field model of the Earth considering J₂ perturbations is of the order of 10⁻⁴ m/s² in terms of acceleration; for low-orbit short-arc segment prediction, it can meet the requirements. The potential function of the earth’s gravity is approximately shown in Equation (2).

$$\mathrm{V}={\displaystyle \frac{\mathrm{GM}}{\mathrm{R}}}\left[1-{\mathrm{J}}_2{\left({\displaystyle \frac{{\mathrm{a}}_{\mathrm{e}}}{\mathrm{R}}}\right)}^2{\mathrm{P}}_2\left(\sin \mathsf{\phi }\right)\right]$$

(2)

In Equation (2), J₂ is the coefficient with a harmonic term, $J_2=0.00108263$, R is the distance from the satellite to the earth’s center, G and M are the earth’s gravitational constant and the earth’s mass, respectively, $G\ast M=3.986005\times {10}^{14} {\mathrm{m}}^3/{\mathrm{s}}^2$, ${\mathrm{a}}_{\mathrm{e}}$ is the Earth’s radius, and φ is the satellite’s geocentric latitude in the earth-solid coordinate system.

$${\mathrm{P}}_2\left(\sin \mathsf{\phi }\right)={\displaystyle \frac{3}{2}}{\left({\displaystyle \frac{\mathrm{z}}{\mathrm{R}}}\right)}^2-{\displaystyle \frac{1}{2}}$$

(3)

In Equation (3), (x,y,z) are the coordinates of the artificial satellite in the earth-solid coordinate system and $\mathrm{R}=\sqrt{{\mathrm{x}}^2+{\mathrm{y}}^2+{\mathrm{z}}^2}$, so the potential function of the earth’s gravitational field is converted into the formula shown in Equation (4).

$$\mathrm{V}={\displaystyle \frac{\mathrm{GM}}{\mathrm{R}}}\left[1-{\mathrm{J}}_2{\left({\displaystyle \frac{{\mathrm{a}}_{\mathrm{e}}}{\mathrm{R}}}\right)}^2\left({\displaystyle \frac{3}{2}}{\left({\displaystyle \frac{\mathrm{z}}{\mathrm{R}}}\right)}^2-{\displaystyle \frac{1}{2}}\right)\right]$$

(4)

By calculating the partial derivative of Equation (4) with respect to x, y, and z, respectively, it can be obtained that the force of the artificial satellite in the earth-solid coordinate system is as shown in Equation (5).

$$F_{\mathrm{Earth}}=\left[\begin{array}{c}\mathrm{GM}\left[-{\displaystyle \frac{1}{{\mathrm{R}}^3}}+{\mathrm{J}}_2{\mathrm{a}}_{\mathrm{e}}^2\left({\displaystyle \frac{15}{2}}{\displaystyle \frac{{\mathrm{z}}^2}{{\mathrm{R}}^7}}-{\displaystyle \frac{3}{2{\mathrm{R}}^5}}\right)\right]\mathrm{x}\\ \mathrm{GM}\left[-{\displaystyle \frac{1}{{\mathrm{R}}^3}}+{\mathrm{J}}_2{\mathrm{a}}_{\mathrm{e}}^2\left({\displaystyle \frac{15}{2}}{\displaystyle \frac{{\mathrm{z}}^2}{{\mathrm{R}}^7}}-{\displaystyle \frac{3}{2{\mathrm{R}}^5}}\right)\right]\mathrm{y}\\ \mathrm{GM}\left[-{\displaystyle \frac{1}{{\mathrm{R}}^3}}+{\mathrm{J}}_2{\mathrm{a}}_{\mathrm{e}}^2\left({\displaystyle \frac{15}{2}}{\displaystyle \frac{{\mathrm{z}}^2}{{\mathrm{R}}^7}}-{\displaystyle \frac{9}{2{\mathrm{R}}^5}}\right)\right]\mathrm{z}\end{array}\right]$$

(5)

Since the position of the artificial satellite is a function of time t, in the earth-solid coordinate system, the position of t at some time near the initial value can be expanded by using Taylor series near the initial value R₀, as shown in Formula (6).

$$R\left(\mathrm{t}\right)=R_0+{R_0}^{\prime }\Delta \mathrm{t}+{\displaystyle \frac{1}{2}}{R_0}^{″}{\Delta \mathrm{t}}^2+{\displaystyle \frac{1}{3!}}{R_0}^{\left(3\right)}{\Delta \mathrm{t}}^3+\dots +{\displaystyle \frac{1}{\mathrm{k}!}}{R_0}^{\left(\mathrm{k}\right)}{\Delta \mathrm{t}}^{\mathrm{k}}+\dots$$

(6)

In Formula (5), take

$$\begin{gathered} \mathrm{A}=\mathrm{GM}\left[-{\displaystyle \frac{1}{{\mathrm{R}}^3}}+{\mathrm{J}}_2{\mathrm{a}}_{\mathrm{e}}^2\left({\displaystyle \frac{15}{2}}{\displaystyle \frac{{\mathrm{z}}^2}{{\mathrm{R}}^7}}-{\displaystyle \frac{3}{2{\mathrm{R}}^5}}\right)\right] \\ \mathrm{B}=\mathrm{GM}\left[-{\displaystyle \frac{1}{{\mathrm{R}}^3}}+{\mathrm{J}}_2{\mathrm{a}}_{\mathrm{e}}^2\left({\displaystyle \frac{15}{2}}{\displaystyle \frac{{\mathrm{z}}^2}{{\mathrm{R}}^7}}-{\displaystyle \frac{3}{2{\mathrm{R}}^5}}\right)\right] \\ \mathrm{C}=\mathrm{GM}\left[-{\displaystyle \frac{1}{{\mathrm{R}}^3}}+{\mathrm{J}}_2{\mathrm{a}}_{\mathrm{e}}^2\left({\displaystyle \frac{15}{2}}{\displaystyle \frac{{\mathrm{z}}^2}{{\mathrm{R}}^7}}-{\displaystyle \frac{9}{2{\mathrm{R}}^5}}\right)\right] \end{gathered}$$

According to Newton’s second law, in the inertial coordinate system of a J2000.0 mean celestial coordinate system, there are formulas as shown in Equation (7).

$$T\ddot{R}=TQR=T\left[\begin{array}{ccc}\mathrm{A}& 0& 0\\ 0& \mathrm{B}& 0\\ 0& 0& \mathrm{C}\end{array}\right]\left[\begin{array}{c}\mathrm{x}\\ \mathrm{y}\\ \mathrm{z}\end{array}\right]$$

(7)

where $T$ is the transformation matrix from the solid coordinate system to the flat celestial coordinate system of J2000.0 and $Q=\left[\begin{array}{ccc}A& 0& 0\\ 0& B& 0\\ 0& 0& C\end{array}\right]$ is a function of the time $t$. It can be seen from $\ddot{R}=QR=\left[\begin{array}{ccc}\mathrm{A}& 0& 0\\ 0& \mathrm{B}& 0\\ 0& 0& \mathrm{C}\end{array}\right]\left[\begin{array}{c}\mathrm{x}\\ \mathrm{y}\\ \mathrm{z}\end{array}\right]$ R that each term in Formula (6) can be written in the form of ${\displaystyle \frac{1}{\mathrm{k}!}}{R_0}^{\left(\mathrm{k}\right)}{\Delta \mathrm{t}}^{\mathrm{k}}={\displaystyle \frac{1}{\mathrm{k}!}}{\Delta \mathrm{t}}^{\mathrm{k}}(F_{\mathrm{k}}{R}_0+G_{\mathrm{k}}\dot{R_0})$, where $F_{\mathrm{k}}$ and $G_{\mathrm{k}}$ are the third-order diagonal square matrices. $F_{\mathrm{k}}$ and $G_{\mathrm{k}}$ can be obtained recursively: $F_0=I$, $G_0=0$, and $F_k={\dot{F}}_{(k-1)}+G_{(k-1)} Q$, $G_k=F_{(k-1)}+{\dot{G}}_{(k-1)}$. If we let $F={\displaystyle {\sum }_{k=0}^{\infty }\frac{1}{k!}}\Delta t^k{F}_k$, $G={\displaystyle {\sum }_{k=0}^{\infty }\frac{1}{k!}}\Delta t^k{G}_k$, and take $R=R_0$, we obtain

$$R\left(\mathrm{t}\right)=F\left(R_0,{\dot{R}}_0,\delta \right)R_0+G(R_0,{\dot{R}}_0,\delta ){\dot{R}}_0$$

(8)

By calculating the terms $F_{\mathrm{k}}$ and $G_{\mathrm{k}}$ from 0 to 6 in sequence, we can calculate $F$ and $G$, and then determine the expression of the satellite’s equation of motion.

Let the coordinate of the target be $(R, A, E)$ in the horizon coordinate system, where $R$ is the distance from the target to the station, in which case

$$\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]=\overline{R} \left[\begin{array}{c}\overline{X}\\ \overline{Y}\\ \overline{Z}\end{array}\right]+\left[\begin{array}{c}{X}_0\\ Y_0\\ Z_0\end{array}\right]$$

(9)

where $(\overline{X},\overline{Y},\overline{Z})$ and $(X_0,Y_0,Z_0)$ are the unit component forms of the observation vector from the station to the satellite in the Earth-fixed coordinate system and the coordinates of the station in the Earth-fixed coordinate system, respectively. Vector $(\overline{X},\overline{Y},\overline{Z})$ is a normalized vector, $\overline{X}=\cos E\ast \cos A$, $\overline{Y}=\cos E\ast \sin A$, $\overline{Z}=\sin E$. According to Equations (8) and (9),

$$F{R}_0+G{\dot{R}}_0=\overline{R}\left[\begin{array}{c}\overline{\mathrm{X}}\\ \overline{\mathrm{Y}}\\ \overline{\mathrm{Z}}\end{array}\right]+\left[\begin{array}{c}{\mathrm{X}}_0\\ {\mathrm{Y}}_0\\ {\mathrm{Z}}_0\end{array}\right]$$

(10)

Let the diagonal elements of diagonal matrix $F$ and $G$ be $(F_1,F_2,F_3)$ and $(G_1,G_2,G_3)$, respectively, then there is

$$\left[\begin{array}{ccc}{F}_1& 0& 0\\ 0& F_2& 0\\ 0& 0& F_3\end{array}\right]\left[\begin{array}{c}{x}_0\\ y_0\\ z_0\end{array}\right]+\left[\begin{array}{ccc}{G}_1& 0& 0\\ 0& G_2& 0\\ 0& 0& G_3\end{array}\right]\left[\begin{array}{c}{\dot{x}}_0\\ {\dot{y}}_0\\ {\dot{z}}_0\end{array}\right]=\overline{R}\left[\begin{array}{c}\overline{X}\\ \overline{Y}\\ \overline{Z}\end{array}\right]+\left[\begin{array}{c}{X}_0\\ Y_0\\ Z_0\end{array}\right]\Rightarrow \left\{\begin{array}{c}{F}_1{x}_0+G_1{\dot{x}}_0=\overline{R}\overline{X}+X_0\\ F_2{y}_0+G_2{\dot{y}}_0=\overline{R}\overline{Y}+Y_0\\ F_3{z}_0+G_3{\dot{z}}_0=\overline{R}\overline{Z}+Z_0\end{array}\right.$$

(11)

And the following basic equation can be obtained by eliminating $\overline{R}$.

$$\left\{\begin{array}{c}{F}_1{x}_0\overline{Y}+G_1{\dot{x}}_0\overline{Y}-F_2{y}_0\overline{X}-G_2{\dot{y}}_0\overline{X}=X_0\overline{Y}-Y_0\overline{X}\\ F_2{y}_0\overline{Z}+G_2{\dot{y}}_0\overline{Z}-F_3{z}_0\overline{Y}-G_3{\dot{z}}_0\overline{Y}=Y_0\overline{Z}-Z_0\overline{Y}\end{array}\right.$$

(12)

An overdetermined system of equations regarding the initial orbit $(x_0,y_0,z_0,{\dot{x}}_0,{\dot{y}}_0,{\dot{z}}_0)$ can be obtained by replacing the observed data at several times into Formula (12), and the initial orbit can be obtained by iteratively solving the overdetermined system. Then, the initial value of the target motion equation is substituted into Formula (8) with the initial orbit as the input time t, and the predicted state value of each moment can be obtained by a numerical solution. The predicted angle value of the corresponding longitude axis and latitude axis can be obtained by converting it to the horizontal coordinate system.

The photoelectric measuring equipment is mainly used for measuring the two-dimensional angle information of the device, namely the beam angle and the angle of the weft axis; the two-axis angle covers the entire upper half of the sphere. Using the current photoelectric measuring equipment of the unit, the total tracking measurement of a low rail satellite is completed, the measurement data are recorded, and the data sampling period is 20 ms. The data are typically collected in 300 s.

The improved Laplacian algorithm was implemented and verified on all the measurement data of a satellite that was actually tracked. The measurement data with a duration of 30 s was used to extrapolate the subsequent orbit’s data, and the extrapolated data and actual measurement data were compared in periods with a low elevation angle and a high elevation angle. The curve of the extrapolated deviation is shown in Figure 1 and Figure 2.

In Figure 1 and Figure 2, the horizontal coordinates are the intervals at which the data were sampled (i.e., the image acquisition period), in units of 20 ms. For example, 1000 in the horizontal coordinates represents the moment at 20 s. The ordinate is the deviation between the predicted target’s angle and the actual measured target’s angle, and the unit is angular seconds. The error of longitude axis and error of latitude axis are the curves of the angular deviation in the longitude and latitude of the targets of the horizontal photoelectric theodolite, respectively. The whole coordinate system reflects the deviation curve of the extrapolated position of the target over time. As time goes on, the deviation gradually diverges. This type of diagram is explained below.

As can be seen from Figure 1 and Figure 2, when the satellite operates at a low elevation angle, this algorithm is used for predicting the orbit, and the angular deviation is small and can be maintained for a long period of time. When the satellite runs at a high elevation angle, the angular deviation predicted by this method is large and diverges more quickly. The reason for this result is that, relative to the center of the Earth, the orbit of the low-orbit satellite is a circular orbit with a fixed angular velocity. In this case, the deviation in the target’s extrapolated position will remain unchanged when converted to the angular deviation. However, compared with the surface observation station, because the position of the satellite moves in an irregular curve relative to the observation station, the angular velocity of its motion is not fixed. At low elevation angles, the angular velocity of its motion is slower, and at high elevation angles, the angular velocity of its motion is constantly faster, as shown in Figure 3. For the extrapolated orbital data with the same deviation in the position, when the satellite is at a low elevation relative to the observation station, the deviation in the predicted position of the target will be reflected in the angular deviation of the horizontal photoelectric theodolite. When the satellite is at a high elevation relative to the observation station, the deviation in the predicted orbital position will be amplified by the angular deviation of the horizontal photoelectric theodolite, as shown in the deviation curve in Figure 1 and Figure 2. This shows that when using the improved Laplace algorithm, it is difficult to guarantee high-precision predictions of the orbit over a long time.

Therefore, the main reason for the low accuracy of the extrapolated data is that the tracking and measurement data used are not accurate enough. In the process of tracking a weak target such as a satellite, an unsteady extraction of the target’s miss distance often occurs, resulting in insufficient accuracy in the measurement data and a large amount of noise. Therefore, this study proposes using the Kalman filter to optimize the estimation of the tracking and measurement data, improve the accuracy of these data, and then improve the accuracy of fitting the equation of motion and the extrapolation data of the orbit.

2.2. Kalman Filtering

The improved Laplacian orbit prediction method requires highly accurate measurement data; otherwise, it is difficult to extrapolate high-precision long-term prediction data. In order to improve the accuracy of the measurement data, the Kalman filter algorithm was used for optimal estimation of the measurement data. Large noise in the measurement data was eliminated to improve the accuracy of the measurement data, thereby reducing the length of the measurement data of tracking and improving the accuracy and time length of the predicted data on the orbit.

2.2.1. Judging the Validity of the Measurement Data

The accuracy of optical measurement equipment in measuring and tracking satellites is mainly reflected in the accuracy of the image processing equipment used to extract the target. As a dim target, a satellite is often extracted incorrectly in the process of image processing. There are two kinds of such errors. One is the extraction error beyond the servo controller’s accuracy and the other is the extraction error within the servo controller’s accuracy.

When the extraction error exceeds the control accuracy of the servo control equipment, there will be a large amount of jitter from the equipment; at the same time, there will be great errors in the recorded measurement data, and the direction of the equipment will deviate from the target trajectory. If the next frame for extraction of the target is correct, the equipment will gradually pull the target back to the center of the field of view, which greatly reduces the stability of tracking. Frequent errors in extracting the target image will lead to a reduction in the accuracy and integrity of the measurement data. For example, as shown in Figure 4, the horizontal coordinates are the points of the sampling interval (i.e., the image acquisition period), and the unit is 20 ms. The ordinate is the miss distance of the target during image processing (i.e., the offset relative to the center of the field of view), in units of 0.66 arcseconds, that is, the pixel resolution of the image.

Assuming that the target can be tracked stably in the first eight frames, the miss distance is 0 in theory, and the error of the miss distance of the target during image extraction is 8 in the ninth frame. If it is in the automatic tracking state at this time, the equipment will turn around, and the target will deviate from the center of the field of view. The target will be extracted correctly again in the tenth frame, and the miss distance will become −8, and then the target will be gradually pulled back to the center of the field of view. In this process, the data recorded in the ninth frame are actually incorrect data. If this phenomenon occurs frequently, it will greatly reduce the stability of tracking and the accuracy and integrity of the measurement data.

In this case, the validity of the measurement data was judged, and the measurement data for tracking were compared with the predicted data of the current frame by using the principle of 3$\sigma$ ($\sigma$ is the precision of the servo controller equipment; the constant term can be appropriately adjusted according to the specific situation). If the deviation is within the range of 3$\sigma$, the tracking data are considered to be valid; otherwise, the tracking data are considered to be invalid and the tracking data of the predicted current frame are directly recorded as the optimal estimated value. The specific formula for judgment is as follows:

$$\left\{\begin{array}{c}\left|L_y-L_c\right|\le 3\sigma \\ \left|B_y-B_c\right|\le 3\sigma \end{array}\right.$$

(13)

where $L_y$ and $B_y$ represent the predicted warp and weft angle values of the current frame’s data and $L_c$ and $B_c$ represent the warp and weft angle values of the current frame’s measured tracking data. When the deviation in the two axes is within the range of 3$\sigma$, the measured tracking data are considered to be valid; otherwise, they are considered to be invalid.

By judging the validity of the measured data, we can eliminate the incorrect data for which the error is greater than the servo controller’s accuracy, as shown in Figure 4.

As can be seen from Figure 5, when the miss distance error of the target extracted from the ninth frame image is 8, the tolerance level of the data determines that these data are wrong, and they are eliminated directly. Prediction of the track extrapolates the current frame’s data as the measurement data and records it and, at the same time, predicts and extrapolates the next frame’s data and guides it. The target always stays in the center of the field of view. When the target is correctly extracted again in the tenth frame, the miss distance is extracted as 0 in theory, and the trajectory of the equipment pointing to the satellite will not shake; even if this phenomenon of error in the extraction of the target occurs frequently, it will not affect the accuracy.

It is necessary to use a filtering method to estimate the optimal measurement data, that is, Kalman filtering, for the measurement data that result from the incorrect extraction of the target image within the precision range of the servo controller.

2.2.2. Kalman Filtering of the Measurement Data

The classical Kalman filtering method was used to filter the measured data within the range of the servo controller’s precision. The classical formula of Kalman filtering is as follows:

$$X(k|k-1)=AX(k-1|k-1)+BU(k)$$

(14)

The prediction state update formula, based on the previous state and control input, predicts the current state at the current time.

$$P\left(k|k-1\right)=AP\left(k-1|k-1\right)A^{\prime }+Q$$

(15)

The prediction covariance update formula is as follows, which updates the prediction state covariance matrix considering the influence of system noise; Q is the system noise.

$$X\left(k|k\right)=X\left(k|k-1\right)+Kg\left(k\right)\left(Z\left(k\right)-HX\left(k|k-1\right)\right)$$

(16)

The state update formula updates the current state based on the measured value and predicted value.

$$Kg\left(k\right)={\displaystyle \frac{P\left(k|k-1\right)H^{\prime }}{\left(HP\left(k|k-1\right)H^{\prime }+R\right)}}$$

(17)

The Kalman gain formula calculates the Kalman gain using the predicted covariance and measurement noise covariance, balances the predicted and measured values, and determines the contribution of the measured values to the state estimation; R is the measured noise.

$$P\left(k|k\right)=\left(I-Kg\left(k\right)H\right)P(k|k-1)$$

(18)

The covariance update formula is used to update the covariance matrix of the state estimation considering the influence of measurement noise.

Kalman filtering is an effective statistic estimation method used to handle systems with noise and uncertainty. The state of the system can be recursively estimated through two steps: prediction and update. The prediction step is based on the dynamic model of the system, while the update step uses actual measurement values to correct the predicted values.

The Kalman filter algorithm used in this study is mainly used for the optimal estimation of the current frame’s data, so the algorithm needs to be improved and simplified in the process of use. The specific improvements are as follows.

(1)

$\mathrm{X}\left(k|k-1\right)$ is the predicted value of the current frame’s data, that is, for the predicted value of the current frame data derived from the improved Laplacian orbit prediction algorithm, Formula (14) is modified to the improved Laplacian orbit prediction algorithm.

(2)

$\mathrm{Z}\left(k\right)$ is the measurement of tracking of the equipment for the current frame’s data, and $\mathrm{X}\left(k|k\right)$ is the optimal estimation of the current frame’s data.

(3)

Because the data on the satellite’s orbit in this study include the azimuth angle and pitch angle, and the two axes of the data are filtered separately, which is appropriate for a single model and a single measurement, Equations (15)–(18) were simplified. The simplified results are as follows.

$$\mathrm{P}\left(k|k-1\right)=\mathrm{P}\left(k-1|k-1\right)+\mathrm{Q}$$

(19)

$$\mathrm{X}\left(k|k\right)=\mathrm{X}\left(k|k-1\right)+\mathrm{Kg}\left(k\right)\left(\mathrm{Z}\left(k\right)-\mathrm{X}\left(k|k-1\right)\right)$$

(20)

$$\mathrm{Kg}\left(k\right)={\displaystyle \frac{\mathrm{P}\left(k|k-1\right)}{\left(\mathrm{P}\left(k|k-1\right)+\mathrm{R}\right)}}$$

(21)

$$\mathrm{P}\left(k|k\right)=\left(1-\mathrm{Kg}\left(k\right)\right)\mathrm{P}(k|k-1)$$

(22)

In the filtering process, the value of $Q$ is set according to the improved Laplacian algorithm for predicting the orbit, and the value of R is set according to the servo controller’s accuracy. Kalman recursion is carried out to filter the orbit data of the current frame. In the process of real-time tracking, filtering is carried out frame by frame, and the optimal estimated value of the filtered tracking data of the current frame is recorded. Precise measurement data are recursively used for improved Laplacian prediction of the orbit to improve the accuracy of prediction.

2.2.3. Workflow for Recording Accurate Measurement Data

The accuracy of the measurement data can be improved by using the method based on the judgments of the validity and Kalman filtering of the measurement data.

As shown in Figure 6, the workflow of recording accurate measurement data includes the following steps:

• Receive the measurement data fed back by the servo control system and the target’s miss distance data fed back by the image processing system and synthesize the target’s measurement data.
• Extrapolate the orbit data of the current frame by using the recorded measurement data and the improved Laplace algorithm for predicting the orbit.
• Compare the target’s measurement data with the extrapolated tracking data of the current frame and judge whether the current frame’s measurement data are valid according to the judgment of validity, that is, whether the target captured by image processing is valid. If it is valid, execute Step 4; if it is invalid, execute Step 5.
• Using the extrapolated tracking data of the current frame in Step 2 as the predicted value and the measured data of the current frame as the measured value, filter the tracking data of the current frame by using the Kalman filtering method to obtain the optimal estimated value of the tracking data of the current frame and record the current frame’s data.
• If several consecutive frames are invalid, the tracking data of the current frame extrapolated in Step 2 are directly used as the optimal estimation value to record the current frame’s data. If consecutive frames are invalid (the number depends on the situation), stop recording.

2.3. An Improved Laplace Satellite Tracking Method Based on the Kalman Filter

The Kalman filter algorithm is used to calculate the optimal estimation of the measurement data and improve the accuracy of the measurement data. Combined with the improved Laplacian algorithm for predicting the orbit, an improved Laplacian satellite tracking method based on the Kalman filter is proposed in this study. Instead of the traditional servo control system, automatic closed-loop tracking based on the interpretation of the target’s miss distance in the image is carried out by using the accurate extrapolated tracking data of the guidance equipment.

2.3.1. Transformation of Coordinates

Since the elevation angle of a satellite’s transit is usually too high, horizontal photoelectric measuring equipment is often used, because it can avoid the zenith’s blind area, which is usually the first choice for the tracking and measurement of satellites. The satellite forecasting data calculated by the SGP4 model is the azimuth and pitch information in horizon mode. Therefore, the coordinates need to be transformed to convert the azimuth and pitch data in horizon mode into longitude and latitude data in horizontal mode to guide the horizontal photoelectric measuring equipment. At the same time, the measurement data of the horizontal photoelectric measuring equipment also need to be converted to the horizon type so that the improved Laplacian prediction of the orbit can be carried out.

The specific formula for transformation of the coordinates is shown below.

(a)

Alt-Az to horizontal is shown in Formula (23):

$$\begin{array}{l}\mathrm{B}=\arctan (\cos \left(\mathrm{E}\right)\ast \sin (\mathrm{A})/\sin \left(\mathrm{E}\right))\\ \mathrm{L}=\arcsin (\cos (\mathrm{A})\ast \cos \left(\mathrm{E}\right))\end{array}$$

(23)

where A and E are the azimuth and pitch angle under the horizontal mode and L and B are the horizontal angles of the longitude and latitude axis after rotation.

(b)

Horizontal to Alt-Az is shown in the Formula (24):

$$\begin{array}{l}\mathrm{A}=\arctan \left(\sin \left(\mathrm{B}\right)/\tan \left(\mathrm{L}\right)\right)\\ \mathrm{E}=\arcsin \left(\cos \left(\mathrm{B}\right)\ast \cos \left(\mathrm{L}\right)\right)\end{array}$$

(24)

where L and B are the angles of the longitude and latitude axes under the horizontal mode and A and E are the converted azimuth and pitch angle.

2.3.2. The Whole Process of the Improved Laplacian Satellite Tracking Method Based on the Kalman Filter

The improved Laplace satellite tracking method based on the Kalman filter is a cyclic and iterative process, including acquisition of the measurement data, extrapolation of the current frame’s data, judgment of the validity and filtering, recording the data, fitting the equation of motion, predicting the next frame’s orbital data, data guidance, and acquisition of the measurement data. The specific process is shown in Figure 7.

The whole process is based on data acquisition, judgments of validity, filtering, and fitting the equation of motion. The specific workflow is broken down as follows:

• The initial unfiltered measurement data are collected. After six frames of data have been obtained, the improved Laplacian algorithm is used for fitting the initial equation of motion, and the next frame’s data are extrapolated.
• The extrapolated data of the next frame are used to guide the server.
• The current frame’s measurement data are collected and extrapolated by the equation of motion, and the validity of the measurement data is determined. If they are valid, Step 4 is performed; if not, Step 6 is performed.
• Kalman filtering is performed on the current frame’s data, the optimal estimation is calculated, and the exact data are recorded, as shown in Step 5.
• If fewer than six frames of accurate data are recorded, the original six frames of data plus the actual recorded data are used to refit the equation of motion. If the actual recorded data are greater than or equal to 6, the equation of motion is refitted using only the actual recorded data. The next frame is extrapolated according to the time step and the moment when the last frame’s data were recorded. Return to Step 2.
• If the number of consecutive invalid frames is small, then perform Step 7; if the number of consecutive invalid frames is large, then perform Step 8 (the number of consecutive invalid frames depends on the situation).
• Directly record the current frame’s extrapolated data as accurate data and return to Step 5.
• Stop recording accurate data and return to Step 5.

3. Results

The improved Laplacian satellite tracking method based on the Kalman filter was implemented on horizontal photoelectric equipment, and the equipment was used to carry out tests. The test was divided into two parts. One involved verifying the effectiveness of the method when the optical measuring equipment lost the target frequently in a short time and lost the target for a long time. The second verified the accuracy of the method for long-term prediction of the orbit and its ability to improve the time.

3.1. Experimental Verification of the Methods

A horizontal photoelectric device was used to track and measure a satellite. The frame frequency of sampling was $50 \mathrm{H}\mathrm{z}$, the pixel resolution of the images was $0.6″$, and the servo controller’s accuracy was better than $2″$. The accuracy of a frame (the current frame) predicted by the improved Laplace algorithm for predicting the orbit based on the Kalman filter was better than $0.5″$, so the Kalman filter constant $Q=0.25, R=4$ was set (this can be adjusted according to the actual situation). We chose a clear day and tried to avoid too small an angle between the satellite’s orbit and the sun to prevent CCD imaging. Two cases of the frequent loss of objects in a short period and a long period were tested.

3.1.1. Frequent Loss of the Target in a Short Period of Time

A satellite with an apparent magnitude of 7.8 was selected during the test. Since the target was dark, unstable extraction of the miss distance often occurred in the process of tracking. In this case, the data on the miss distance obtained by the improved Laplacian satellite tracking method based on the Kalman filter are shown in Figure 8 and Figure 9.

As can be seen from the figures, there is often a large jump in miss distance, which is caused by the error of the miss distance in image extraction because the target is dark. If the automatic tracking method of image processing is used at this time, the target may be lost and the miss distance will diverge, and it is difficult for the miss distance to converge again. Even if the predicted number is switched immediately, the miss distance will jitter back and forth for a short time, and then gradually converge back to the center of the field of view. When using this method to guide and track the satellite’s data throughout the whole process, the program will make a judgment according to the predicted data of the current frame and the measured data fed back by the equipment. If the deviation exceeds 3σ, the data will be directly eliminated and the predicted data will be stored as the optimal estimate of the current frame’s data. Meanwhile, the next frame’s data on the orbit will be predicted for guidance. At this time, even if the image processing has a large jump in the miss distance, the device will follow the correct trajectory through guidance, and the target will still be in the center of the field of view. Namely, the improved Laplacian satellite tracking method based on the Kalman filter can significantly improve the stability of the tracking process and ensure the integrity and accuracy of the measurement data.

3.1.2. Losing the Target for Long Periods of Time

During the test, a satellite was selected, and its apparent magnitude was 6.5. When the device accurately tracked the satellite for a period of time (more than 30 s), it manually switched the miss distance to an invalid state on the program interface. At this time, the program would automatically switch the tracking measurement data to the invalid state and stop recording the accurate measurement data. The best estimate of the recorded orbital data is used to predict all subsequent data. Because it is invalid to switch the miss distance manually at this time, the miss distance of the target actually exists all the time. After a period of time, the miss distance is switched to the effective state again, and the target is always near the center of the field of view. After recording the initial measurement data of more than five frames, the program automatically restarts the forecast data guidance tracking and pulls the target into the center of the field of view again to continue the tracking and measurement. The specific test data are shown in Figure 10 and Figure 11.

After the equipment stably tracked the target for a period of time, that is, after the optimal estimate of the satellite orbit data was stably obtained for a period of time, setting the miss distance was invalid, and the complete predictive data-guided tracking was performed. As can be seen from the figure, the deviation of miss distance gradually increases at this time. After about 2000 frames, it is effective to set the miss distance again. At this time, the deviation of miss distance is about 7, and the target is still near the center of the field of view. After a few frames of initial measurement data acquisition, the program automatically restarts the orbit prediction data guidance, and soon the target is pulled back to the center of the field of view to continue the tracking measurement. The method proposed in this paper can predict long-term and high-precision orbit data even if the target is lost for a long time and ensure that the target is near the center of the field of view. When the target appears again after occlusion, accurate tracking and measurement work can be continued.

3.2. Test Verification of Orbit Prediction Capability

The same horizontal photoelectric measuring equipment is used to track and measure the satellite and record all the measurement data. During the test, a clear evening was selected so that CCD could not be imaged due to the small angle with the sun to avoid affecting the collection of test data. The long time prediction of the satellite orbit was carried out without and with the Kalman filter, and the orbit prediction capability of the proposed method was verified. The satellite measured in this experiment is a low-orbit sun-synchronous orbit satellite, and multiple sets of data were recorded. Due to space problems, three sets of data are used below for comparison. The usage of measurement data is shown in Figure 12.

Figure 10 and Figure 11 describe the difference curve between the predicted orbit data and the real measured orbit data. Th e position of horizontal coordinate 0 is the vertex of the target; the measurement data 30 s before the vertex were used for motion model calculation to predict the next data, and the difference was compared with the real measured orbit data. The magnitude and divergence of the deviation value repre-sent the ability of orbit prediction.

3.2.1. When Kalman Filtering Was Not Used

The improved Laplacian algorithm for predicting the orbit was used to forecast the orbits of multiple satellites for a long time without effective judgment and Kalman filtering of the measurement data. The curve of the deviation in the predictions is shown in Figure 13.

The above three tests obtained all the predicted deviation data during the high elevation period of satellite transit. The tracking measurement data 30 s before the vertex were used to fit the motion equation, and then the extrapolated orbit prediction data were used to make a difference with the actual measured data to obtain the deviation curve. According to the data results of the three tests, for the extrapolated data for 20 s, the maximum deviation is about$4″$.

3.2.2. When Kalman Filtering Was Used

The improved Laplacian satellite tracking method based on the Kalman filter proposed in this study was used to carry out long-term prediction of the orbits of multiple satellites under the conditions of effective judgment of the measurement data and Kalman filtering. The curve of deviation in the prediction is shown in Figure 14.

The above three tests also obtained the prediction deviation data during the period of a high elevation angle of satellite transit. The measurement data 30 s before the vertex were used, and after effective data judgment and Kalman filtering, the motion equation was fitted. Then, the deviation curve was obtained by using the extrapolated orbit prediction data and the actual measured data. For the extrapolated data for 30 s, the maximum deviation is less than 3$″$.

3.2.3. Comparison of the Data

Due to the limitations of space, multiple sets of data from the two tests were analyzed, and the statistical data are shown in

Table 1. Statistical table of experimental data from Section 3.2.1.

Number	Length of Data	Length of Extrapolated Data	Maximum Deviation	Number of Counts
1	$20 \mathrm{s}$	$20 \mathrm{s}$	$7.27″$	26
2	$30 \mathrm{s}$	$20 \mathrm{s}$	$4.61″$	22
3	$40 \mathrm{s}$	$20 \mathrm{s}$	$2.75″$	14

and

Table 2. Statistical table of experimental data from Section 3.2.2.

Number	Length of Data	Length of Extrapolated Data	Maximum Deviation	Number of Counts
1	$20 \mathrm{s}$	$30 \mathrm{s}$	$4.66″$	27
2	$25 \mathrm{s}$	$30 \mathrm{s}$	$3.58″$	21
3	$30 \mathrm{s}$	$30 \mathrm{s}$	$2.81″$	16

. Table 1 shows the case without judgment of the validity and Kalman filtering of the measured data using the improved Laplacian algorithm for predicting the orbit. Table 2 shows the data obtained by the improved Laplacian satellite tracking method based on Kalman filtering proposed in this study.

The statistics in Table 1 show that without the judgment of validity and Kalman filtering of the measurement data, the length of the extrapolated high-precision tracking data was generally in the order of 20 s, and measurement data with a length of more than 40 s were required to ensure that the maximum deviation of the predicted data was less than $3^{″}$.

As can be seen from the statistical data in Table 2, the length of the high-precision extrapolated tracking data could reach 30 s by using the optimal estimates of the tracking data recorded after filtering to fit the equation of motion and predict the track. If the maximum deviation of the predicted data was less than $3″$, the length of the used measuring data was less than 30 s.

It can be seen from the analysis results of the two groups of test data that the duration of the predicted data is significantly increased and the duration of the measured data is significantly reduced by using the Kalman filter compared with not using the Kalman filter. With a prediction deviation of less than $3″$, the extrapolation time was increased by 50% from 20 s to 30 s, and the length of the measurement data used was reduced by 25% from 40 s to 30 s.

4. Conclusions

In this study, an improved Laplacian satellite tracking method based on the Kalman filter is proposed. The method performs judgments of the validity and Kalman filtering of the measured data in the process of real-time tracking and measurement and calculates the optimal estimate of the current frame’s orbital data, which improves the accuracy of the measured orbital data. The best estimate of the measured tracking data is used to improve the Laplacian method of fitting the orbit’s equation of motion, and the data on the orbit are extrapolated frame by frame to guide the equipment to realize the whole process of tracking the satellite.

At present, this method is used for a certain type of equipment. An analysis of the test data showed that the method did not deviate from the satellite’s trajectory, and the target was still in the center of the field of view when the target was lost frequently in a short time (i.e., the miss distance of the target could not be extracted), which ensured the stability of the tracking equipment. For the case when the target could not be extracted for a long time, this method was used to predict data to guide the device. The target was stable near the center of the field of view for a long time. When the target appeared again, (i.e., when the miss distance of the target could be stably extracted again), the tracking measurement was carried out again, realizing blind tracking in a real sense. Moreover, the test data showed that judgment of the validity and filtering of the measured data significantly improved the length of time and accuracy of the extrapolated tracking data and reduced the length of the measured tracking data required. Under the condition that the maximum deviation of the predicted data is better than $3″$, the length of the extrapolated data increased from the original $20 \mathrm{s}$ to $30 \mathrm{s}$. The length of the measured tracking data required reduced from $40 \mathrm{s}$ to $30 \mathrm{s}$. The conclusion is that the improved Laplacian satellite tracking method based on the Kalman filter can obviously improve the stability of the tracking process and ensure the integrity and accuracy of the measured data.