Short arc initial orbit determining method based on gauss solution cluster

Abstract

The invention discloses a preferable short arc initial orbit determining method based on Gauss solution cluster, and belongs to the astrodynamics field, including: grouping the observation data, using Gauss method to obtain the target state vector at the corresponding time point for each group of data, forming a solution set of preliminary estimation; dividing the solution set of preliminary estimation into a position component vector solution set and a velocity component vector solution set for clustering to obtain a position component vector solution cluster and a velocity component vector solution cluster; based on the position component vector solution cluster and the velocity component vector solution cluster, generating a two-dimensional trajectory solution set; evaluating each of the two-dimensional trajectories by using a trajectory optimal method, calculating the number of root of orbits corresponding to the optimal two-dimensional trajectory, thereby completing determination of initial orbit.

Claims

1. A short arc initial orbit determining method for single space-based imaging observation platform based on Gauss solution cluster, wherein the method comprises the following steps: S1. dividing an observation data, which are captured by an observation camera, into groups, wherein for each group of data, Gauss method is adopted to find a target state vector for a corresponding time point to form a solution set of preliminary estimation; S2. dividing a solution set of preliminary estimation into a position component vector solution set and a velocity component vector solution set, which are respectively clustered to obtain a position component vector solution cluster and a velocity component vector solution cluster; S3. generating a two-dimensional trajectory solution set based on the position component vector solution cluster and the velocity component vector solution cluster; S4. using a trajectory optimization method to evaluate each two-dimensional trajectory, calculate the number of roots of orbit corresponding to the optimal two-dimensional trajectory, and complete determination of the short arc initial orbit; wherein step S2 comprises the following sub-steps: S201. eliminating non-reasonable solutions from the solution set of preliminary estimation; S202. dividing the remaining state vector solution set into a position component vector solution set and a velocity component vector solution set, and clustering being performed respectively so correct solutions are gathered in the same cluster as many as possible, thereby obtaining a position component vector solution cluster and a velocity component vector solution cluster; S203. performing noise reduction preprocessing on the position component vector solution cluster and the velocity component vector solution cluster, respectively.

2. The short arc initial orbit determining method according to claim 1, wherein step S1 comprises the following sub-steps: S101. dividing the observation data into groups, and vector data corresponding to every 3 time points being sorted into a group; S102. for each group of vector data, using the conventional Gauss method to find a target state vector at the corresponding time point to form a solution set of preliminary estimation.

3. The short arc initial orbit determining method according to claim 2, wherein it is assumed that an observation time point is [t.sub.0,t.sub.k], and k+1≥3, and for input data [t.sub.0,t.sub.k] of a total of k+1 time points, the interval parameter Nrate=└rate*k┘ is taken, wherein └•┘ represents the rounding down, rate is the interval rate, and the grouped observation data is expressed as {{right arrow over (L)}.sub.t.sub.s1, {right arrow over (L)}.sub.t.sub.s2, {right arrow over (L)}.sub.t.sub.s3, {right arrow over (R)}.sub.t.sub.s1, {right arrow over (R)}.sub.t.sub.s2, {right arrow over (R)}.sub.t.sub.s3}, wherein {right arrow over (R)}.sub.t.sub.s represents a satellite observation position vector, {right arrow over (L)}.sub.t.sub.s represents an observation angle vector $.Math.\frac{k-\mathrm{Nrate}}{2}.Math.\le s2\le .Math.\frac{k-\mathrm{Nrate}}{2}.Math.,s1=.Math.\frac{2*s2+\mathrm{Nrate}-k+1}{2}.Math.,s3=.Math.\frac{2*s2+\mathrm{Nrate}-k+1}{2}.Math..$

4. The short arc initial orbit determining method according to claim 1, wherein position data and velocity data satisfying any one of conditions |{right arrow over (r)}.sub.t.sub.s|<6378.14 km or |{right arrow over (r)}.sub.t.sub.s|>42378 km or |{right arrow over (v)}.sub.t.sub.s|>11.2 km/s are excluded, wherein {right arrow over (r)}.sub.t.sub.s represents the target position vector at time point t.sub.s, and {right arrow over (v)}.sub.t.sub.s represents the target velocity vector at the time point t.sub.s.

5. The short arc initial orbit determining method according to claim 1, wherein in step S202, a k-means clustering method is used for clustering and Chauvenet's—criterion discriminating method is adopted to eliminate abnormal data.

6. The short arc initial orbit determining method according to claim 1, wherein step S3 comprises the following sub-steps: S301. performing fitting on the position component vector solution cluster and the velocity component vector solution cluster after noise reduction to construct a state vector combination at various time points; S302. generating a three-dimensional trajectory solution set of the target orbit according to the state vector combination corresponding to various time points; S303. projecting the 3D trajectory solution set of the target orbit onto an instantaneous observation image plane according to a measurement status of an observation platform to obtain a 2D trajectory solution set.

7. The short arc initial orbit determining method according to claim 1, wherein step S4 comprises the following sub-steps: S401. calculating a derivative error and a position error for each two-dimensional trajectory; wherein the formula for calculating the derivative error of the m-th two-dimensional trajectory is as follows: $\Delta {\mathrm{uv}}^{\prime m}=\underset{n=0}{\overset{k}{.Math.}}\left[.Math.\left(2t_n\times {\hat{U}}^{″m}+{\hat{U}}^{\prime m}\right)-\left(2t_n\times U^{″*}+U^{\prime *}\right).Math.+.Math.\left(2t_n\times {\hat{V}}^{″m}+{\hat{V}}^{\prime m}\right)-\left(2t_n\times V^{″*}+V^{\prime *}\right).Math.\right]{\hat{U}}_{t_n}^m={\hat{U}}^{″m}{t}_n^2+{\hat{U}}^{\prime m}{t}_n+{\hat{U}}_0^m{\hat{V}}_{t_n}^m={\hat{V}}^{″m}{t}_n^2+{\hat{V}}^{\prime m}{t}_n+{\hat{V}}_0^m{U}_{t_n}^{*}=U^{″*}{t}_n^2+U^{\prime *}{t}_n+U_0^{*}{V}_{t_n}^{*}=V^{″*}{t}_N^2+V^{\prime *}{t}_n+V_0^{*}$
.Math..sub.t.sub.n.sup.m=.Math.″.sup.mt.sub.n.sup.2+.Math.′.sup.mt.sub.n+.Math..sub.0.sup.m
{circumflex over (V)}.sub.t.sub.n.sup.m={circumflex over (V)}″.sup.mt.sub.n.sup.2+{circumflex over (V)}′.sup.mt.sub.n+{circumflex over (V)}.sub.0.sup.m
U.sub.t.sub.n*=U″*t.sub.n.sup.2+U′*t.sub.n+U.sub.0*
v.sub.t.sub.n*=V″*t.sub.n.sup.2+V′*t.sub.n+V.sub.0* wherein U″*, V″*, U′*, V′* are quadratic polynomial fitting parameters of the two-dimensional trajectory {l*=(U.sub.t.sub.n, V.sub.t.sub.n*, t.sub.n, n=0, 1, 2, 3 . . . k)}, .Math..sub.t.sub.n.sup.m,{circumflex over (V)}.sub.t.sub.n.sup.m is a corresponding image plane coordinates of an estimated trajectory l.sub.m′ at a time point t.sub.n, m represents the m-th two-dimensional trajectory, and k is the total number of observation time points; the formula for calculating the position error of the m-th two-dimensional trajectory is as follows: ${\mathrm{Rt}}^m=\underset{n=0}{\overset{k}{.Math.}}\frac{D_{t_n}^m}{{\mathrm{Dc}}_{t_n}^m}{D}_{t_n}^m=\sqrt{{\left({\hat{U}}_{t_n}^m-U_{t_n}^{*}\right)}^2+{\left({\hat{V}}_{t_n}^m-V_{t_n}^{*}\right)}^2}{\mathrm{Dc}}_{t_n}^m=\sqrt{{\left(U_c^m-U_{t_n}^{*}\right)}^2+{\left(V_c^m-V_{t_n}^{*}\right)}^2}$
D.sub.t.sub.n.sup.m=√{square root over ((.Math..sub.t.sub.n.sup.m−U.sub.t.sub.n*).sup.2+({circumflex over (V)}.sub.t.sub.n.sup.m−V.sub.t.sub.s*).sup.2)}
Dc.sub.t.sub.n.sup.m=√{square root over ((U.sub.c.sup.m−U.sub.t.sub.n*).sup.2+(V.sub.c.sup.m−V.sub.t.sub.n*).sup.2)} wherein, D.sub.t.sub.n.sup.m is an image plane distance between a predicted value and an observed value, Dc.sub.t.sub.n.sup.m is the image plane distance between the intersection point of the estimated trajectory and the observed trajectory and the observed point, t.sub.n is the observation time, .Math..sub.t.sub.n.sup.m,{circumflex over (V)}.sub.t.sub.n.sup.m is the corresponding image plane coordinate of the estimated trajectory l.sub.m′ at the time point t.sub.n, U.sub.t.sub.n*V.sub.t.sub.n* is the corresponding image plane coordinate of the observed trajectory l* at the time point t.sub.n, and U.sub.c.sup.mV.sub.c.sup.m is the image plane coordinate of the intersection point of the estimated trajectory l.sub.m′ and the observed trajectory; S402. taking the derivative error and the position error into comprehensive consideration to select the optimal two-dimensional trajectory, and the number of root of orbit corresponding thereto is used as an optimal estimated value; S403. outputting the number of root of orbit corresponding to the optimal two-dimensional trajectory, and the determination of the short arc initial orbit is completed.

8. The short arc initial orbit determining method according to claim 7, wherein step S402 is specifically as follows: for the m-th orbit, the derivative error and position error are sorted from small to large, respectively, and the ranking sequence numbers Rank.sub.vel,m and Rank.sub.dis,m are obtained, the two sequence numbers are between 1˜N.sub.sv, and the sum of the two is calculated as the final error evaluation:
Rank.sub.m=Rank.sub.vel,m+Rank.sub.dis,m for the final error evaluation, the smallest one m.sub.opt is taken as the optimal predicted trajectory. $m_{\mathrm{opt}}=\underset{m}{\arg \min }\left\{{\mathrm{Rank}}_m;m=1,2,3.Math.N_{\mathrm{sv}}\right\}$ wherein, N.sub.sv is the number of orbit root obtained by the combination of the target position and the velocity data corresponding to various time points.

9. A computer-readable storage medium, wherein the computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, the short arc initial orbit determining method claimed in claim 1 is implemented.

10. A computer-readable storage medium, wherein the computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, the short arc initial orbit determining method claimed in claim 2 is implemented.

11. A computer-readable storage medium, wherein the computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, the short arc initial orbit determining method claimed in claim 3 is implemented.

12. A computer-readable storage medium, wherein the computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, the short arc initial orbit determining method claimed in claim 4 is implemented.

13. A computer-readable storage medium, wherein the computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, the short arc initial orbit determining method claimed in claim 5 is implemented.

14. A computer-readable storage medium, wherein the computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, the short arc initial orbit determining method claimed in claim 6 is implemented.

15. A computer-readable storage medium, wherein the computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, the short arc initial orbit determining method claimed in claim 7 is implemented.

16. A computer-readable storage medium, wherein the computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, the short arc initial orbit determining method claimed in claim 8 is implemented.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) FIG. 1 is a flowchart of a short arc initial orbit determining method based on a Gauss solution cluster according to an embodiment of the present disclosure.

(2) FIG. 2 is a schematic diagram showing space-based observation according to an embodiment of the present disclosure.

(3) FIG. 3 is a solution cluster diagram of a target position obtained by the Gauss method according to an embodiment of the present disclosure.

(4) FIG. 4 is a clustering diagram of k-means clustering performed on position solution cluster according to an embodiment of the present disclosure.

(5) FIG. 5 is a position solution clustering diagram after the adopted Savitzky-Golay is smoothed according to an embodiment of the present disclosure.

(6) FIG. 6 is a three-dimensional diagram corresponding to various position solutions obtained by calculating orbit roots according to an embodiment of the present disclosure.

(7) FIG. 7 is a two-dimensional diagram showing projection of the three-dimensional orbit corresponding to various position solutions projected on the image plane according to an embodiment of the present disclosure.

(8) FIG. 8 is a three-dimensional diagram of an optimal orbit and an actual observation orbit according to an embodiment of the present disclosure.

(9) FIG. 9 is a flowchart of sub-steps of step S2 based on a Gauss solution cluster according to an embodiment of the present disclosure.

(10) FIG. 10 is a flowchart of sub-steps of step S3 based on a Gauss solution cluster according to an embodiment of the present disclosure.

(11) FIG. 11 is a flowchart of sub-steps of step S4 based on a Gauss solution cluster according to an embodiment of the present disclosure.

DESCRIPTION OF THE EMBODIMENTS

(12) In order to make the objectives, technical solutions, and advantages of the present disclosure clearer, the following further describes the present disclosure in detail with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present disclosure and are not intended to limit the present disclosure.

(13) First of all, some terms related to the present disclosure are explained as follows.

(14) Solution Cluster: A Collection of Solutions

(15) As shown in FIG. 1, the present disclosure provides a preferred short arc initial orbit determining method based on a Gauss solution cluster. The method includes the following steps:

(16) S1. The observation data is divided into groups, and for each group of data, the Gauss method is adopted to find the target state vector for the corresponding time point to form a solution set of preliminary estimation.

(17) S2. The solution set of preliminary estimation is divided into a position component vector solution set and a velocity component vector solution set, which are respectively clustered to obtain a position component vector solution cluster and a velocity component vector solution cluster.

(18) S3. A two-dimensional trajectory solution set is generated based on the position component vector solution cluster and the velocity component vector solution cluster.

(19) S4. The trajectory optimization method is adopted to evaluate each two-dimensional trajectory, calculate the number of root of orbit corresponding to the optimal two-dimensional trajectory, and complete the determination of the initial orbit.

(20) Step S1. The observation data is divided into groups, and for each group of data, the Gauss method is adopted to find the target state vector for the corresponding time point to form a solution set of preliminary estimation.

(21) S101. The observation data is divided into groups, and the vector data corresponding to every 3 time points is sorted into a group.

(22) The observation data is observation data from space-based monocular imaging platform or Earth observation data. The observation data includes: the observation position (the position of the platform itself relative to the geocentric coordinate system) corresponding to each observation time point, the observation angle vector, satellite observation position vector represented by variant {right arrow over (R)}.sub.t.sub.s at time point t.sub.s, and observation angle vector represented by {right arrow over (T)}.sub.t.sub.s at time point t.sub.s. The position of the observation satellite, the angle of the observation camera, and the observation imaging data at the same time point constitute a pair of data, and multiple time point data constitute a group of data.

(23) In order to construct a solution cluster subsequently, the present disclosure first needs to group the observation data, and the vector data corresponding to every 3 time points is sorted into a group. Assume that the observation time point is [t.sub.0,t.sub.k], and k+1≥3, and for input data [t.sub.0,t.sub.k] of a total of k+1 time points, the interval parameter Nrate=└rate*k┘ is taken, wherein └•┘ represents the rounding down, rate is the interval rate, normally the taken value is 0.2, and the grouped observation data is expressed as {{right arrow over (L)}.sub.t.sub.s1, {right arrow over (L)}.sub.t.sub.s2, {right arrow over (L)}.sub.t.sub.s3, {right arrow over (R)}.sub.t.sub.s1, {right arrow over (R)}.sub.t.sub.s2, {right arrow over (R)}.sub.t.sub.s3}, wherein {right arrow over (R)}.sub.t.sub.s represents a satellite observation position vector, {right arrow over (L)}.sub.t.sub.s represents an observation angle vector,

(24) $.Math.\frac{k-\mathrm{Nrate}}{2}.Math.\le s2\le .Math.\frac{k-\mathrm{Nrate}}{2}.Math.,s1=.Math.\frac{2*s2+\mathrm{Nrate}-k+1}{2}.Math.,s3=.Math.\frac{2*s2+\mathrm{Nrate}-k+1}{2}.Math..$

(25) S102. For each group of vector data, the conventional Gauss method is adopted to find the target state vector at the corresponding time point to form a solution set of preliminary estimation.

(26) Since the conventional Gauss method can only input three data and the input information is limited, the present disclosure pondered on using the Gauss method for multiple times to increase data utilization. The number of root of short arc target orbit finally determined based on a combination of multiple sets of data will be closer to the true value.

(27) The Gauss method solves the target position vector and the target velocity vector at the intermediate time point s2 by inputting the observation position vector and the observation angle vector corresponding to three different time points. The target refers to an unpowered target, which can be a debris or a runaway satellite, etc. According to the method of the present disclosure, every three time points in a series of time points are sorted into a group, and each group of data is solved by the conventional Gauss method once to obtain the target position vector and velocity vector solution at the intermediate time point in the group of data. The number of solutions obtained from each group of vector data may be 1, 2, and 3. The Gauss method is adopted multiple times to find the target vector corresponding to different time points to form a collection of solutions of preliminary estimation.

(28) Step S2. The solution set of preliminary estimation is divided into a position component vector solution set and a velocity component vector solution set, which are respectively clustered to obtain a position component vector solution cluster and a velocity component vector solution cluster.

(29) S201. The non-reasonable solutions in the solution set of preliminary estimation are eliminated.

(30) The position data and velocity data satisfying any one of conditions |{right arrow over (r)}.sub.t.sub.s|<6378.14 km (radius of Earth) or |{right arrow over (r)}.sub.t.sub.s|>42378 km (geostationary orbit radius) or |{right arrow over (v)}.sub.t.sub.s|>11.2 km/s (second cosmic velocity) are excluded, wherein {right arrow over (r)}.sub.t.sub.s represents the target position vector at time point t.sub.s, and {right arrow over (v)}.sub.t.sub.s represents the target velocity vector at the time point t.sub.s.

(31) Obviously, a solution that satisfies any of the above three conditions is a non-reasonable solution. Therefore, the target position data in which the distance between the to-be-solved target position and center of Earth is smaller than the radius of the Earth or larger than the geostationary orbit radius and the target velocity data greater than the second cosmic velocity are excluded.

(32) S202. The remaining state vector solution set is divided into a position component vector solution set and a velocity component vector solution set, and clustering is performed respectively so that correct solutions are gathered in the same cluster as many as possible, thereby obtaining the position component vector solution cluster and the velocity component vector solution cluster.

(33) In order to separate the correct root from the wrong root so that the correct roots are gathered in the same cluster as many as possible, the position component vector solution set and the velocity component vector solution set are clustered respectively. The k-means clustering method is preferably adopted by the present disclosure for clustering (k is a value of 2 or 3), and Chauvenet's—criterion is adopted for discrimination to exclude abnormal data. The abnormal data refers to data that exhibits abnormality in the same group, such as outliers.

(34) If the number of position solutions corresponding to various time points is less than or equal to 1, k-means clustering is not required; if the maximum number of position solutions corresponding to any time point is 2, then k-means clustering where k=2 is performed once; if the maximum number of position solutions corresponding to any time point is 3, a k-means clustering where k=3 is performed once. The velocity component vector solution set is clustered in the same way as described above.

(35) S203. Noise reduction preprocessing is performed on the position component vector solution cluster and the velocity component vector solution cluster, respectively.

(36) The present disclosure performs noise reduction preferably based on a Savitzky-Golay filtering method.

(37) After all roots are processed by removing non-reasonable data and smoothing noise reduction, the accuracy of trajectory optimization can be improved.

(38) Step S3. A two-dimensional trajectory solution set is generated based on the position component vector solution cluster and the velocity component vector solution cluster.

(39) S301. The position component vector solution cluster and the velocity component vector solution cluster after the noise reduction are fitted to construct a state vector combination at various time points.

(40) The data after noise deduction is utilized to perform fitting on the target position component vector and the target velocity component vector. The combination of velocity vector and position vector corresponding to the time point t.sub.s is ({right arrow over (r)}.sub.t.sub.s′, {right arrow over (v)}.sub.t.sub.s′, t.sub.s).

(41) S302. A three-dimensional trajectory solution set of the target orbit is generated according to the state vector combination corresponding to various time points.

(42) With the combination ({right arrow over (r)}.sub.t.sub.s′, {right arrow over (v)}.sub.t.sub.s′, t.sub.s) of the velocity vector and position vector corresponding to a certain time point t.sub.s, the number of root {Ô.sub.m=(a.sub.m, e.sub.m, i.sub.m, Ω.sub.m, ω.sub.m, M.sub.m)} of orbit corresponding to ({right arrow over (r)}.sub.t.sub.s′, {right arrow over (v)}.sub.t.sub.s′, t.sub.s) is calculated. Based on the number of root of orbit, it can be determined that the three-dimensional trajectory solution set is {{circumflex over (l)}.sub.m=({circumflex over (x)}.sub.t.sub.n.sup.m,ŷ.sub.t.sub.n.sup.m,{circumflex over (z)}.sub.t.sub.n.sup.m, t.sub.n)}, wherein a is the semi-major axis, e is the first eccentricity, i is the inclination of the orbit, Q is the right ascension of the ascending intersection, ω is the perigee angle, and M is the mean anomaly, ({circumflex over (x)}.sub.t.sub.n.sup.m,{circumflex over (x)}.sub.t.sub.n.sup.m,{circumflex over (x)}.sub.t.sub.n.sup.m) is the corresponding three-dimensional coordinate of the observed trajectory {circumflex over (l)}.sub.m at the time point t.sub.n, m=1, 2, . . . , N.sub.sv, n=0, 1, . . . , k, N.sub.sv is the total number of root of orbit obtained based on the combination of the target positions and the velocity data corresponding to various time points.

(43) S303. The 3D trajectory solution set of the target orbit is projected to the instantaneous observation image plane according to the measurement status of the observation platform to obtain a 2D trajectory solution set.

(44) The 3D trajectory solution set {{circumflex over (l)}.sub.m=({circumflex over (x)}.sub.t.sub.n.sup.m,ŷ.sub.t.sub.n.sup.m,{circumflex over (z)}.sub.t.sub.n.sup.m,t.sub.n)} of the target orbit is projected to the instantaneous observation image plane, thereby obtaining the 2D trajectory solution set {l.sub.m′=.Math..sub.t.sub.n.sup.m,{circumflex over (V)}.sub.t.sub.n.sup.m, t.sub.n}, wherein {circumflex over (l)}.sub.m is m-th 3D observed trajectories, ({circumflex over (x)}.sub.t.sub.n.sup.m,{circumflex over (x)}.sub.t.sub.n.sup.m,{circumflex over (x)}.sub.t.sub.n.sup.m) is the corresponding 3D coordinate of the observed trajectory {circumflex over (l)}.sub.m at the time point t.sub.n, l.sub.m′ is the m-th estimated trajectory, .Math..sub.t.sub.n.sup.m,{circumflex over (V)}.sub.t.sub.n.sup.m is the corresponding image plane coordinate of the estimated trajectory l.sub.m′ at the time point t.sub.n, n=0, 1, . . . , k, m=1, 2, . . . , N.sub.sv.

(45) Step S4. The trajectory optimization method is adopted to evaluate each two-dimensional trajectory, calculate the number of roots of orbit corresponding to the optimal two-dimensional trajectory, and complete the determination of the initial orbit.

(46) S401. The derivative error and position error for each two-dimensional trajectory are calculated.

(47) The formula for calculating the derivative error of the m-th two-dimensional trajectory is as follows:

(48) $\Delta {\mathrm{uv}}^{m\prime }=\underset{n=0}{\overset{k}{.Math.}}\left[.Math.\left(2t_n\times {\hat{U}}^{″m}+{\hat{U}}^{\prime m}\right)-\left(2t_n\times U^{″*}+U^{\prime *}\right).Math.+.Math.\left(2t_n\times {\hat{V}}^{″m}+{\hat{V}}^{\prime m}\right)-\left(2t_n\times V^{″*}+V^{\prime *}\right).Math.\right]{\hat{U}}_{t_n}^m={\hat{U}}^{″m}{t}_n^2+{\hat{U}}^{\prime m}{t}_n+{\hat{U}}_0^m{\hat{V}}_{t_n}^m={\hat{V}}^{″m}{t}_n^2+{\hat{V}}^{\prime m}{t}_n+{\hat{V}}_0^m{U}_{t_n}^{*}=U^{″*}{t}_n^2+U^{\prime *}{t}_n+U_0^{*}{V}_{t_n}^{*}=V^{″*}{t}_n^2+V^{\prime *}{t}_n+V_0^{*}$

(49) wherein U″*, V″*, U′*, V′* are the quadratic polynomial fitting parameters of the two-dimensional trajectory {l*=(U.sub.t.sub.n, V.sub.t.sub.n*, t.sub.n, n=0, 1, 2, 3 . . . k)}, .Math..sub.t.sub.n.sup.m,{circumflex over (V)}.sub.t.sub.n.sup.m is the corresponding image plane coordinates of the estimated trajectory l.sub.m′ at the time point t.sub.n, m represents the m-th two-dimensional trajectory, and k is the total number of observation time points.

(50) The camera parameters (orientation) remain the same, and the phase plane coordinates (U.sub.t.sub.n*, V.sub.t.sub.n*) of various time points with respect to the observation of target can be obtained. Quadratic polynomial fitting is performed on the observation data for two times in U and V directions respectively, thereby obtaining the second term and the fourth term.

(51) The formula for calculating the position error of the m-th two-dimensional trajectory is as follows:

(52) ${\mathrm{Rt}}^m=\underset{n=0}{\overset{k}{.Math.}}\frac{D_{t_n}^m}{{\mathrm{Dc}}_{t_n}^m}{D}_{t_n}^m=\sqrt{{\left({\hat{U}}_{t_n}^m-U_{t_n}^{*}\right)}^2+{\left({\hat{V}}_{t_n}^m-V_{t_n}^{*}\right)}^2}{\mathrm{Dc}}_{t_n}^m=\sqrt{{\left(U_c^m-U_{t_n}^{*}\right)}^2+{\left(V_c^m-V_{t_n}^{*}\right)}^2}$

(53) wherein, D.sub.t.sub.n.sup.m is the image plane distance between the predicted value and the observed value, Dc.sub.t.sub.n.sup.m is the image plane distance between the intersection point of the estimated trajectory and the observed trajectory and the observed point, t.sub.n is the observation time, .Math..sub.t.sub.n.sup.m,{circumflex over (V)}.sub.t.sub.n.sup.m is the corresponding image plane coordinate of the estimated trajectory l.sub.m′ at the time point t.sub.n, U.sub.t.sub.n*V.sub.t.sub.n* is the corresponding image plane coordinate of the observed trajectory l* at the time point t.sub.n, and U.sub.c.sup.mV.sub.c.sup.m is the image plane coordinate of the intersection point of the estimated trajectory l.sub.m′ and the observed trajectory.

(54) S402. The derivative error and the position error are taken into comprehensive consideration to select the optimal two-dimensional trajectory, and the number of root of orbit corresponding thereto is used as the optimal estimated value.

(55) For the m-th orbit, the derivative error and position error are sorted from small to large, respectively, and the ranking sequence numbers Rank.sub.vel,m and Rank.sub.dis,m are obtained, the two sequence numbers are between 1˜N.sub.sv, and the sum of the two is calculated as the final error evaluation.
Rank.sub.m=Rank.sub.vel,m+Rank.sub.dis,m

(56) For the final error evaluation, take the smallest value m.sub.opt as the optimal two-dimensional trajectory:

(57) $m_{\mathrm{opt}}=\underset{m}{\arg \min }\left\{{\mathrm{Rank}}_m;m=1,2,3.Math.N_{\mathrm{sv}}\right\}$

(58) wherein, N.sub.sv is the number of root of orbits obtained based on the combination of the target position and the velocity data corresponding to various time points.

(59) S403. The number of root of orbit corresponding to the optimal two-dimensional trajectory is output, and the determination of the initial orbit is completed.

(60) The present disclosure is suitable for single space-based imaging observation platform, long-range unpowered targets under short-term observation conditions (greater than 5 s), such as debris, out-of-control satellites, etc. The disclosure can effectively solve the problem of multiple roots that cannot be solved by the conventional Gauss method, and has higher utilization of original observation data. Moreover, the present disclosure adopts the concept of solution cluster and optimal method to improve the accuracy of the initial orbit determination. Under the same input error conditions, the technical solution of the present disclosure has better performance than other schemes in terms of the short-term observation data, that is, the present disclosure is more suitable for short-arc observation than other disclosures.

Example 1

(61) As shown in Table 1, Table 1 shows data of observation satellites and number of root of target orbits in Example 1.

(62) TABLE-US-00001 TABLE 1 a (km) e i (°) u (°) Ω (°) Observing 6978.1 0 97.6893 0.1719 275.5927 space- based data Target 6987.0 7.8e−09 79.1255 4.7555 272.7279 data

(63) S101. The observation data is divided into groups, and the vector data corresponding to every 3 time points is sorted into a group.

(64) With respect to the input data [t.sub.0, t.sub.170] of a total of 171 time points, the interval parameter Nrate=└rate*k┘ is taken, the observation data {{right arrow over (L)}.sub.t.sub.s1, {right arrow over (L)}.sub.t.sub.s2, {right arrow over (L)}.sub.t.sub.s3, {right arrow over (R)}.sub.t.sub.s1, {right arrow over (R)}.sub.t.sub.s2, {right arrow over (R)}.sub.t.sub.s3} is sorted into a group, wherein └•┘ represents the rounding down, rate is the interval rate, and

(65) $\mathrm{rate}=0.2,s1=.Math.\frac{2*s2+34-170+1}{2}.Math.,.Math.\frac{170-34}{2}.Math.\le s2\le .Math.\frac{170+34}{2}.Math.,s3=.Math.\frac{2*s2-34+170-1}{2}.Math.,$
{right arrow over (R)}.sub.t.sub.s represents a satellite observation position vector, {right arrow over (L)}.sub.t.sub.s represents an observation angle vector.

(66) With respect to the observation satellites and the number of root of target orbits in the data in Table 1, simulation is conducted to obtain the data regarding observation satellites and partial position of the target as shown in Table 2.

(67) TABLE-US-00002 TABLE 2 Point {right arrow over (r)}.sub.t (km) {circumflex over (L)}.sub.t (incl. error/km) {right arrow over (R)}.sub.t (km) number t (s) x.sub.t.sub.s y.sub.t.sub.s z.sub.t.sub.s x.sub.t.sub.s y.sub.t.sub.s z.sub.t.sub.s x.sub.t.sub.s y.sub.t.sub.s z.sub.t.sub.s 0 0.0 306.30 −6980.14 −52.03 −0.9877 −0.0934 −0.1254 705.00 −6942.43 −1.41 . . . . . . . . . 34 3.4 311.15 −6980.06 −26.81 −0.9872 −0.09439 −0.1285 701.58 −6942.73 24.01 . . . . . . . . . 68 6.8 315.99 −6979.90 −1.59 −0.9867 −0.0954 −0.1318 698.15 −6942.94 49.44 . . . . . . . . . 102 10.2 320.84 −6979.64 23.62 −0.9861 −0.0964 −0.1351 694.71 −6943.06 74.87 . . . . . . . . . 136 13.6 325.67 −6979.28 48.84 −0.9855 −0.0976 −0.1387 691.26 −6943.09 100.30 . . . . . . . . . 170 17.0 330.37 −6978.85 73.31 −0.9848 −0.0987 −0.1424 687.80 −6943.01 125.73

(68) S102. For each group of vector data, the conventional Gauss method is adopted to find the target state vector at the corresponding time point to form a solution set of preliminary estimation.

(69) A single group of data after grouping is as shown in FIG. 2, point E represents the center of the coordinate system. The distance between the observation point and the target is defined as ρ.sub.i, then:

(70) 0$\{\begin{array}{c}\stackrel{.fwdarw.}{r_1}=\stackrel{.fwdarw.}{R_1}+{\rho }_1\stackrel{.fwdarw.}{L_1}\\ \stackrel{.fwdarw.}{r_2}=\stackrel{.fwdarw.}{R_2}+{\rho }_2\stackrel{.fwdarw.}{L_2}\\ \stackrel{.fwdarw.}{r_3}=\stackrel{.fwdarw.}{R_3}+{\rho }_3\stackrel{.fwdarw.}{L_3}\end{array}$

(71) wherein, the satellite observation position vector at the time point t.sub.i and the corresponding position vector of the target are {right arrow over (R)}.sub.i and {right arrow over (r)}.sub.i respectively, and the unit sight vector corresponding to the observation time points t.sub.1, t.sub.2 and t.sub.3 are {right arrow over (L)}.sub.1, {right arrow over (L)}.sub.2 and {right arrow over (L)}.sub.3 respectively. There are six unknown vectors in the above formula, namely, {right arrow over (r)}.sub.1, {right arrow over (r)}.sub.2 and {right arrow over (r)}.sub.3 as well as ρ.sub.1, ρ.sub.2 and ρ.sub.3; r.sub.2 is the magnitude of {right arrow over (r)}.sub.2, and the obtained target position is as shown in FIG. 3, the target velocity is as shown in FIG. 4, wherein the mark * indicates the calculated position solution, the mark + indicates the actual position solution of target, and the mark .box-tangle-solidup. indicates the observation satellite position. According to the Gauss method, the initial value of {right arrow over (r)}.sub.2 and {right arrow over (v)}.sub.2 can be solved, and the data obtained from this example is shown in Table 3. Table 3 shows data concerning target position and velocity solution cluster based on Gauss method in Example 1.

(72) TABLE-US-00003 TABLE 3 position vector velocity vector Position vector data Combination Point.sub.t.sub.s1 solution set based solution set based on after removal of non- of Point.sub.t.sub.s2 on Gauss method Gauss method (m/s, reasonable data (km, sequence Point.sub.t.sub.s3 (km, km, km) m/s, m/s) km, km)  1 Point 0, (−162380.97, (1037557.63, 47310.82, (256.55, −6985.65, Point 68, −22712.67, −19460.97); −9.53); Point −21732.39); (1801.66, 79.89, (880.40, −6925.32, 136 (256.55, 7407.1059); 73.78) −6985.65, −9.53); (−2171.29, −101.28, (880.40, 7510.17) −6925.32, 73.78) . . . . . . . . . . . . . . . 34 Point (−145427.83, (948881.44, 43250.12, (247.73, −6986.78, 34, −21233.37, −17276.82); 12.88); Point −19943.24); (1892.90, 110.84, (874.68, −6925.47, 102, (247.73, −6986.78, 7404.08); 98.77) Point 12.88); (−2182.56, 170 (874.68, −6925.47, −74.81, 7510.29) 98.77)

(73) S201. The non-reasonable solutions in the solution set of preliminary estimation are eliminated.

(74) The principle for elimination is that positions and velocity data that satisfy |{right arrow over (r)}.sub.t.sub.s|<6378.14 km, |{right arrow over (r)}.sub.t.sub.s|>42378 km and |{right arrow over (v)}.sub.t.sub.s|>11.2 km/s are excluded. The obtained positions and velocity in the example are as shown in Table 4, which shows the target position error data before smoothing in Example 1.

(75) TABLE-US-00004 TABLE 4 Optimized Simulated Position Point.sub.t.sub.s1 position position vector Point.sub.t.sub.s2 vector vector error Set Point.sub.t.sub.s3 (km, km, km) (km, km, km) (km, km, km)  1 Point 0, (256.55, (315.99, (59.44, 5.74, Point 68, −6985.65, −6979.90, 7.94) Point 136 −9.53) −1.59) (−564.40, −54.57, (880.40, −75.38) −6925.32, 73.78) . . . . . . . . . . . . . . . 34 Point 34, (247.73, (320.84, (73.10, 7.14, Point 102, −6986.78, −6979.64, 10.73) Point 170 12.88) 23.62) (−553.84, −54.16, (874.68, −75.14) −6925.47, 98.77)

(76) It can be expressed as follows.

(77) $\left\{\left({\stackrel{.fwdarw.}{r}}_{t_s}^m,t_s\right);m=1,2,3.Math.N;\left(.Math.\frac{170-34}{2}.Math.\le s\le .Math.\frac{170+34}{2}.Math.\right)\right\}\left\{\left({\stackrel{.fwdarw.}{v}}_{t_s}^m,t_s\right);m=1,2,3.Math.N;\left(.Math.\frac{170-34}{2}.Math.\le s\le .Math.\frac{170+34}{2}.Math.\right)\right\}$

(78) wherein, N is the number of remaining positions or velocity solutions after removing non-reasonable data.

(79) S202. The remaining state vector solution set is divided into a position component vector solution set and a velocity component vector solution set, and clustering is performed respectively so that correct solutions are gathered in the same cluster as many as possible, thereby obtaining the position component vector solution cluster and the velocity component vector solution cluster.

(80) Perform k-means clustering on the obtained

(81) $\left\{\left({\stackrel{.fwdarw.}{r}}_{t_s}^m,t_s\right),\left({\stackrel{.fwdarw.}{v}}_{t_s}^m,t_s\right);m=1,2,3.Math.N;\left(.Math.\frac{170-34}{2}.Math.\le s\le .Math.\frac{170+34}{2}.Math.\right)\right\}$
respectively. Preferably, if the number of position solutions corresponding to all the time points is 1, k-means clustering is not required; if the maximum number of position solutions corresponding to any time point is 2, then k-means clustering where k=2 is performed once; if the maximum number of position solutions corresponding to any time point is 3, a k-means clustering where k=3 is performed once. The velocity component vector solution set is clustered in the same way as described above.

(82) Taking the location solution as an example, the clustering result is shown in FIG. 5, wherein the mark * indicates cluster 1, the mark o indicates cluster 2, the mark + is the target actual positions, and the mark .box-tangle-solidup. indicates the observation satellite position. The data in this example is as shown in Table 5, which shows the target position error data in Example 1 after smoothing.

(83) TABLE-US-00005 TABLE 5 .sup.a (km) .sup.e incl (°) .sup.ω (°) Ω (°) M (°) .sup.u (°) Target 6987.0 7.8e−09 79.1255 / 272.7279 / 4.7555 data Data of 6928.366 0.008 80.6249 185.7025 272.7467 173.7632 −0.5343 present method Original −0.236 143393.4 84.1812 88.1617 2.7930 0 88.1617 solution data 1 based on Gauss method Original 8447.614 0.1748 112.7304 11.4240 278.8715 352.7453 4.1693 solution data 2 based on Gauss method Original 6805.225 0.0262 86.0354 189.0760 273.7521 170.6514 −0.2726 solution data 3 based on Gauss method

(84) Example 1 adopts Chauvenet's —criterion for discrimination to eliminate abnormal data.

(85) When the confidence probability is taken as

(86) $1-\frac{1}{2n},$
the confidence probability formula of Chauvenet's—criterion can be expressed as: w.sub.n=1+0.4In(n).

(87) With respect to the value x.sub.i to be tested, the absolute value of the difference between the value x.sub.i and the sample mean is calculated. If the following formula is satisfied, the current value to be tested is excluded: |x.sub.i−x|>w.sub.nS.sub.x, wherein x.sub.i is the value to be tested. Here, the components in the three directions (x, y, z) of the position and velocity vectors are selected for three detections respectively. If any one of the detections satisfies the above formula, then it is excluded. x is the average value of such sample, w.sub.n is the confidence probability of the Chauvenet's criterion, S.sub.x is the standard deviation of such sample.

(88) S203. Noise reduction preprocessing is performed on the position component vector solution cluster and the velocity component vector solution cluster, respectively.

(89) The present disclosure performs noise reduction preferably based on a Savitzky-Golay filtering method, and the processing is as follows:

(90) ${\left(\begin{array}{ccccc}{\stackrel{\_}{x}}_{t_{s-2}}& {\stackrel{\_}{x}}_{t_{s-1}}& {\stackrel{\_}{x}}_{t_s}& {\stackrel{\_}{x}}_{t_{s+1}}& {\stackrel{\_}{x}}_{t_{s+1}}\\ {\stackrel{\_}{y}}_{t_{s-1}}& {\stackrel{\_}{y}}_{t_{s-1}}& {\stackrel{\_}{y}}_{t_s}& {\stackrel{\_}{y}}_{t_{s+1}}& {\stackrel{\_}{y}}_{t_{s+1}}\\ {\stackrel{\_}{z}}_{t_{s-1}}& {\stackrel{\_}{z}}_{t_{s-1}}& {\stackrel{\_}{z}}_{t_s}& {\stackrel{\_}{z}}_{t_{s+1}}& {\stackrel{\_}{z}}_{t_{s+1}}\end{array}\right)}^T=\frac{1}{70}\left(\begin{array}{ccccc}69& 4& -6& 4& -1\\ 4& 54& 24& -16& 4\\ -6& 24& 34& 24& -6\\ 4& -16& 24& 54& 4\\ -1& 4& -6& 4& 69\end{array}\right){\left(\begin{array}{ccccc}{x}_{t_{s-2}}& x_{t_{s-1}}& x_{t_s}& x_{t_{s+1}}& x_{t_{s+2}}\\ y_{t_{s-2}}& y_{t_{s-1}}& y_{t_s}& y_{t_{s+1}}& y_{t_{s+2}}\\ z_{t_{s-2}}& z_{t_{s-1}}& z_{t_s}& z_{t_{s+1}}& z_{t_{s+2}}\end{array}\right)}^T$

(91) After using the above-mentioned five-point and three-order template for smoothing, the noise can be effectively reduced in three directions (x, y, z). The position data after smoothing will be more concentrated in the solution space. Take the position solution cluster as an example, the results before and after smoothing are shown in FIG. 6, wherein the mark * indicates the target position before smoothing, the mark o indicates the target position after smoothing.

(92) S301. The position component vector solution cluster and the velocity component vector solution cluster after the noise reduction are fitted to construct a state vector combination at various time points.

(93) For the data after smoothing and noise-reduction, a quadratic polynomial fitting is performed in three directions of x, y, and z. After calculating the estimated values of the target positions corresponding to each intermediate time point

(94) $\left\{t_s^{\prime };.Math.\frac{170-34}{2}.Math.\le s\le .Math.\frac{170+34}{2}.Math.\right\}$
on the x, y, and z components, they can be combined to obtain the corresponding position correction estimated value

(95) $\left\{{{\stackrel{.fwdarw.}{r}}_{t_s}}^{\prime };.Math.\frac{170-34}{2}.Math.\le s\le .Math.\frac{170+34}{2}.Math.\right\}$
and velocity correction estimated value

(96) $\left\{{{\stackrel{.fwdarw.}{r}}_{t_s}}^{\prime };.Math.\frac{170-34}{2}.Math.\le s\le .Math.\frac{170+34}{2}.Math.\right\}.$

(97) S302. A three-dimensional trajectory solution set of the target orbit is generated according to the state vector combination corresponding to various time points.

(98) S303. The 3D trajectory solution set of the target orbit is projected to the instantaneous observation image plane according to the measurement status of the observation platform to obtain a 2D trajectory solution set.

(99) Find the corresponding target fitting positions at the same time points and combine them to obtain 2×2×34 pairs of positions and velocity and time combination ({right arrow over (r)}.sub.t.sub.s′, {right arrow over (v)}.sub.t.sub.s′, t.sub.s). Calculate the corresponding number of root of orbits {Ô.sub.m=(a.sub.m, e.sub.m, incl.sub.m, Ω.sub.m, ω.sub.m, M.sub.m, Ω.sub.m); m=1, 2, 3, . . . , N.sub.sv}, the corresponding target orbits {{circumflex over (l)}.sub.m=({circumflex over (x)}.sub.t.sub.n.sup.m,ŷ.sub.t.sub.n.sup.m,{circumflex over (z)}.sub.t.sub.n.sup.m, t.sub.n, n=0, 1, 2, 3 . . . k); m=1, 2, 3, . . . , N.sub.sv}, and the trajectory of the target orbit projected onto the image plane: {l*=(U.sub.t.sub.n, V.sub.t.sub.n*, t.sub.n, n=0, 1, 2, 3 . . . k); m=1, 2, 3, . . . , N.sub.sv}. The obtained target predicted orbit is as shown in FIG. 7, wherein the mark indicates an estimated orbit, and mark o indicates a target simulated orbit.

(100) S401. The derivative error and position error for each two-dimensional trajectory are calculated.

(101) The formula for calculating the derivative error of the m-th two-dimensional trajectory is as follows:

(102) $\Delta {\mathrm{uv}}^{\prime m}=\underset{n=0}{\overset{k}{.Math.}}\left[.Math.\left(2t_n\times {\hat{U}}^{″m}+{\hat{U}}^{\prime m}\right)-\left(2t_n\times U^{″*}+U^{\prime *}\right).Math.+.Math.\left(2t_n\times {\hat{V}}^{″m}+{\hat{V}}^{\prime m}\right)-\left(2t_n\times V^{″*}+V^{\prime *}\right).Math.\right]$${\hat{U}}_{t_n}^m={\hat{U}}^{″m}{t}_n^2+{\hat{U}}^{\prime m}{t}_n+{\hat{U}}_0^m{\hat{V}}_{t_n}^m={\hat{V}}^{″m}{t}_n^2+{\hat{V}}^{\prime m}{t}_n+{\hat{V}}_0^m{U}_{t_n}^{*}=U^{″*}{t}_n^2+U^{\prime *}{t}_n+U_0^{*}{V}_{t_n}^{*}=V^{″*}{t}_n^2+V^{\prime *}{t}_n+V_0^{*}$

(103) wherein U″*, V″*, U′*, V′* are the quadratic polynomial fitting parameters of the two-dimensional trajectory {l*=(U.sub.t.sub.n, V.sub.t.sub.n*, t.sub.n, n=0, 1, 2, 3 . . . k)}, .Math..sub.t.sub.n.sup.m,{circumflex over (V)}.sub.t.sub.n.sup.m is the corresponding image plane coordinates of the estimated trajectory l.sub.m′ at the time point t.sub.n, m represents the m-th two-dimensional trajectory, and k is the total number of observation time points.

(104) The camera parameters (orientation) remain the same, and the phase plane coordinates (U.sub.t.sub.n*, V.sub.t.sub.n*) of various time points with respect to the observation of target can be obtained. Quadratic polynomial fitting is performed on the observation data for two times in U and V directions respectively, thereby obtaining the second term and the fourth term.

(105) The formula for calculating the position error of the m-th two-dimensional trajectory is as follows:

(106) ${\mathrm{Rt}}^m=\underset{n=0}{\overset{k}{.Math.}}\frac{D_{t_n}^m}{{\mathrm{Dc}}_{t_n}^m}{D}_{t_n}^m=\sqrt{{\left({\hat{U}}_{t_n}^m-U_{t_n}^{*}\right)}^2+{\left({\hat{V}}_{t_n}^m-V_{t_n}^{*}\right)}^2}{\mathrm{Dc}}_{t_n}^m=\sqrt{{\left(U_c^m-U_{t_n}^{*}\right)}^2+{\left(V_c^m-V_{t_n}^{*}\right)}^2}$

(107) wherein, D.sub.t.sub.n.sup.m is the image plane distance between the predicted value and the observed value, Dc.sub.t.sub.n.sup.m is the image plane distance between the intersection point of the estimated trajectory and the observed trajectory and the observed point, t.sub.n is the observation time, .Math..sub.t.sub.n.sup.m,{circumflex over (V)}.sub.t.sub.n.sup.m is the corresponding image plane coordinate of the estimated trajectory l.sub.m′ at the time point t.sub.n, U.sub.t.sub.n*V.sub.t.sub.n* is the corresponding image plane coordinate of the observed trajectory l* at the time point t.sub.n, and U.sub.c.sup.mV.sub.c.sup.m is the image plane coordinate of the intersection point of the estimated trajectory l.sub.m′ and the observed trajectory.

(108) S402. The derivative error and the position error are taken into comprehensive consideration to select the optimal two-dimensional trajectory, and the number of root of orbit corresponding thereto is used as the optimal estimated value.

(109) For the m-th orbit, the derivative error and position error are sorted from small to large, respectively, and the ranking sequence numbers Rank.sub.vel,m and Rank.sub.dis,m are obtained, the two sequence numbers are between 1˜N.sub.sv, and the sum of the two is calculated as the final error evaluation.
Rank.sub.m=Rank.sub.vel,m+Rank.sub.dis,m

(110) For the final error evaluation, take the smallest value m.sub.opt as the optimal two-dimensional trajectory:

(111) 0$m_{\mathrm{opt}}=\underset{m}{\arg \min }\left\{{\mathrm{Rank}}_m;m=1,2,3.Math.N_{\mathrm{sv}}\right\}$

(112) wherein, N.sub.sv is the number of root of orbits obtained based on the combination of the target position and the velocity data corresponding to various time points.

(113) S403. The number of root of orbit corresponding to the optimal two-dimensional trajectory is output, and the determination of the initial orbit is completed.

(114) As shown in FIG. 8, the mark represents all estimated orbits, .box-tangle-solidup. represents a preferred estimated orbit, and mark o is a target simulated orbit.

(115) The above are only the preferred embodiments of this application, but the scope of this application is not limited thereto. Any changes or replacements that can be easily derived by person skilled in the art based on the technical scope disclosed in the present disclosure should fall within the technical scope disclosed in this application. Therefore, the scope sought to be protected by this application shall be subject to the scope of the claims.