SEARCH AND TRACKING METHOD FOR FULL TIME-DOMAIN LASER DETECTION OF SPACE DEBRIS, STORAGE MEDIUM, AND ELECTRONIC DEVICE

Abstract

In a search and tracking method for full time-domain laser detection of space debris, a set of latest precision orbital parameters of a debris object and start and end moments of a current transit of the object are first obtained. Search-specific guidance data is generated based on the above information and in combination with estimation of a maximum along-track error of the orbital parameters of the object during the current transit. A DLR system performs multi-elevation search on the object based on the search-specific guidance data, obtains a plurality of pieces of detection data of the object after detecting the object during the search, determines an along-track error of the orbital parameters of the object based on the detection data, and corrects the orbital parameters of the object in real time based on the along-track error, so as to guide the DLR system to subsequently track and detect the object.

Description

TECHNICAL FIELD

The present disclosure relates to the field of space exploration, and specifically, to a search and tracking method for full time-domain laser detection of space debris, a storage medium, and an electronic device.

BACKGROUND

At the beginning of this century, to ensure the safety of an in-orbit spacecraft during operation, increasing attention is paid to collision avoidance and active clearance technologies for space debris (including defunct satellites, rocket wreckage, and debris generated by object disintegration, collision, explosion, and the like). Effective implementation of these technologies depends on accurate acquisition and prediction of the orbit of the space debris. However, the accuracy of debris orbit determination based on radar and optical measurement data are difficult to meet the requirements of those technologies. Therefore, a debris laser ranging (DLR) technology has become a new development direction.

Essentially, there is no difference between the DLR technology and a satellite laser ranging (SLR) technology in terms of the detection principle, and both the DLR technology and the SLR technology obtain the range of a space object by measuring the round-trip time of a laser signal from a ground station to the space object. However, the DLR technology faces a greater technical challenge and difficulty due to different detection objects, and there are two key technical problems that need to be resolved. The first problem is related to detection of echo signals. The detection object of the SLR technology is a cooperative object. For a laser beam emitted by a system, the vast majority of photons can be reflected back by retroreflectors of the object according to the original path, making it easy to detect the echo signals. However, the detection object of the DLR technology is space debris (non-cooperative object). For the laser beam emitted by the system, only a very small proportion of photons can return to the ground station by diffuse reflection from target surface, making it difficult to detect the echo signals. The second problem is related to generation of echo signals, and the laser beam emitted by the system is required to be able to detect (hit) the debris and generate echo signals, thereupon. This is usually achieved through observation guidance provided by orbit prediction of the object. Because the laser beam is extremely narrow, high accuracy is required for the orbit prediction. For a small quantity of cooperative space objects targeted by the SLR technology, a high-precision ephemeris file (CPF) updated and published daily by the International Laser Ranging Service (ILRS) can be used to achieve sufficient accurate orbital predictions of the object to guide a SLR system to effectively track and measure the object. However, for a general debris object targeted by the DLR technology, such a high-precision ephemeris is unavailable.

In order to make the laser beam to detect and track the debris object, the optical guidance method is currently used worldwide. Firstly, an optical telescope or astronomical camera with a large field of view is used to search and detect the object to obtain a certain amount of angle observation data. Then the laser beam is guided to detect and track the object (see non-patent documents 1-3).

In consideration of the dependence of current debris laser detection on an optical device, application of the DLR technology is greatly limited for the following reasons. First, optical measurement is generally constrained by skylight and Earth shadow conditions, which limits the DLR technology to only being effectively implemented around twilight (within these few hours, space debris on orbit is illuminated by the sun while the ground laser ranging station is still in darkness). In addition, due to being limited by laser energy, at present, the DLR technology can only effectively detect the low earth orbit (LEO) debris with a range of no more than 3,000 kilometers. Since it usually takes a long time for the low earth orbit (LEO) debris to pass through the Earth's shadow, after daylight at the laser station and the time the object passing through the Earth's shadow are deducted, it is estimated that a total time window available for DLR detection does not exceed 6 hours (see non-patent document 4). This not only causes a great waste of detection resources, but also results in unsatisfactory orbital phase distribution of measurement data, which is not conducive to effective improvement of the accuracy of debris orbit determination, thus seriously reducing application efficacy of DLR technology. In order to effectively expand the time window of DLR detection, a daylight space debris laser detection technology has been explored internationally in recent years (see non-patent document 4). Implementation of the daylight space debris laser detection technology mainly benefits from the improvement of optical observation technology. The daylight space debris laser detection technology still relies on optical guidance and does not fully eliminate the impact of noise, such as a skylight background. Therefore, only a small quantity of close-range debris with large size could be detected by daylight laser ranging currently, which means the application efficacy of DLR technology is still significantly limited. With further improvement of the optical observation technology, even if the capability of space debris daylight laser detection is significantly improved in the future, the time window problem within the DLR detection still exists due to optical invisibility of a debris object in the Earth's shadow region.

In summary, the time window problem has become a main technical bottleneck that restricts the application efficacy of the current DLR technology. The time window is generated due to initial guidance of optical observation during implementation of the DLR technology. However, for the SLR technology, due to the guidance of high-precision orbital predictions, there is no time window problem in the application process, and laser detection can be achieved 24 hours a day. Therefore, the generation of the time window is not necessarily related to the DLR technology itself. If the debris object can still be detected and tracked without relying on guidance of the optical device, the DLR system will fundamentally break through the limitation of the time window, and its application efficacy will be essentially improved. In view of this, the present disclosure proposes a search and tracking method that can achieve full time-domain laser detection of the space debris, along with the corresponding storage medium and electronic device. Referring to the Chinese patent with a publication number “CN113640787A” and titled “Constant Elevation Search Method for Space Object Capturing Using Narrow Beam Radar” (see patent document 1), a detection process of the debris search and tracking method is still divided into a constant elevation search process and a tracking process after detecting the object. However, since laser beam for detection is much narrower than radar's (a beam diameter of the DLR system generally does not exceed 20 arc seconds, while a beam diameter of the narrow beam radar is usually not less than 200 arc seconds, in other words, there is a difference of magnitude between the beam diameter of the DLR system and the beam diameter of the narrow beam radar), the DLR system and the narrow beam radar have significantly different technical requirements for object search and tracking.

CITED DOCUMENTS

Non-Patent Documents

• • 1. Bin Li, Jizhang Sang and Zhongping Zhang, A Real-Time Orbit Determination Method for Smooth Transition from Optical Tracking to Laser Ranging of Debris, Sensors, 16, 962, 2016.
  • 2. Michael A. Steindorfer, Georg Kirchner, Franz Koidl, Peiyuan Wang, Alfredo Anton, Jaime Fernandez, Diego Escobar, Jiri Silha, Klaus Merz, Stare and Chase: Optical Pointing Determination, Orbit Calculation and Satellite Laser Ranging Within a Single Pass, Proc. 7th European Conference on Space Debris, Darmstadt, Germany, 18-21 Apr. 2017.
  • 3. J. Rodriguez, E. Cordelli, and T. Schildknecht, The new Stare & Chase Procedure at The Swiss Optical Ground Station and Geodynamics Observatory Zimmerwald, Proc. 1st NEO and Debris Detection Conference, Darmstadt, Germany, 22-24 Jan. 2019.
  • 4. Michael A. Steindorfer, Georg Kirchner, Franz Koidl, Peiyuan Wang, Beatriz Jilete & Tim Flohrer, Daylight space debris laser ranging, Nature communications, 2020.

Patent Documents

• • 1. Chinese patent with a publication number “CN113640787A” and titled “Constant Elevation Search Method for Space Object Capturing Using Narrow Beam Radar” on 12 Nov. 2021.

SUMMARY

In order to overcome the limitations in the prior art, the present disclosure provides a search and tracking method for full time-domain laser detection of space debris, a storage medium, and an electronic device, with main objectives of fully expanding effective detection resources of a DLR system and optimizing a phase distribution structure of DLR data, thereby greatly improving orbit determination and prediction accuracy of the space debris. Further, the improvement of orbit determination and prediction accuracy provides powerful technical support for implementation of collision avoidance and active clearance tasks of the space debris, thereby significantly improving a safety assurance capability of an in-orbit operating spacecraft.

To achieve the above objectives, the present disclosure adopts following technical solutions:

A search and tracking method for full time-domain laser detection of space debris is implemented during each transit of a space debris object relative to a DLR system, and includes:

• • based on a set of precision orbital parameters of an object that are determined based on data of precision radars and/or charge-coupled device (CCD) telescopes of a space monitoring network, and by taking an intermediate moment determined based on pre-obtained start and end moments of a current transit of the object as an initial time, generating a set of initial quasi-mean elements, and constructing an analytical perturbation model;
  • discretizing an along-track error, obtaining, in combination with the analytical perturbation model, a series of new perturbation models for describing motion of all virtual objects, and forming an apparent path cluster by different apparent paths generated from the virtual objects during the current transit;
  • determining a maximum detectable elevation h_q and a minimum detectable elevation h_p of each valid apparent path in the apparent path cluster, and then calculating search elevations of the DLR system;
  • for each valid apparent path in the apparent path cluster, calculating a detection elevation based on a corresponding perturbation model and in consideration of an impact of the light travel time, and obtaining, by using a dichotomous approach based on a difference between the detection elevation and a search elevation, detection feature parameters of virtual objects corresponding to all apparent paths when the virtual objects rise to the search elevation;
  • generating a series of sets of search-specific guidance data based on the detection feature parameters of the virtual objects, where each set of search-specific guidance data corresponds to one dwell of a laser beam emitted by the DLR system on the search elevation;
  • performing, by the DLR system, multi-elevation search on the object based on the search-specific guidance data; and
  • switching to, by the DLR system, a tracking process after detecting the object, calculating an along-track error of the initial quasi-mean elements based on a few sets of detection data obtained before switching to the tracking process, to correct the along-track error of the initial quasi-mean elements, and generating tracking-specific guidance data based on a set of corrected initial quasi-mean elements to guide the DLR system to track the object.

In order to optimize the above technical solution, following specific measures are also used.

Further, after constructing the analytical perturbation model, the method further includes: calculating, in combination with the constructed analytical perturbation model, a moment of a theoretical closest approach point of the object during the current transit; and for the apparent path cluster formed by the different apparent paths generated from the virtual objects during the current transit, determining the visibility of each apparent path in the apparent path cluster based on the calculated moment of the theoretical closest approach point, removing all the invisible apparent paths and further removing all the apparent paths with an unsatisfactory detectable arc length, and taking the remaining apparent paths as valid apparent paths; and

• • after obtaining the detection feature parameters of the virtual objects corresponding to all the valid apparent paths when the virtual objects rise to the search elevation, performing continuation processing on a value of an azimuth in the detection feature parameters to make the value of the azimuth change continuously.

Further, the constructing an analytical perturbation model specifically includes:

• • selecting the intermediate moment T₀ of the current transit of the object as a reference moment, defining an orbital along-track error of the object during the current transit as an along-track error at the reference moment, and using the first-type nonsingular orbital elements as basic variables to construct an analytical perturbation model with the reference moment as the initial time:
  • assuming that one set of known precision orbital parameters of the object are t_q, {right arrow over (r)}_q, {right arrow over ({dot over (r)})}_q, and ϵ_a, where t_q represents the epoch moment of the one set of orbital parameters, {right arrow over (r)}_q and {right arrow over ({dot over (r)})}_q respectively represent the position vector and the velocity vector of the object relative to an epoch geocentric inertial coordinate system, and ϵ_a represents area-mass ratio of the object, performing perturbation propagation by using a numerical method and a precision dynamic model based on the known orbital parameters t_q, {right arrow over (r)}_q, {right arrow over ({dot over (r)})}_q, and ϵ_a, performing from the moment t_q to the moment T₀ to obtain a position vector {right arrow over (r)}₀ and a velocity vector {right arrow over ({dot over (r)})}₀ of the object relative to the epoch geocentric inertial coordinate system at the moment T₀, performing conversion on {right arrow over (r)}₀ and {right arrow over ({dot over (r)})}₀ to obtain a set of initial quasi-mean elements of the object at the moment T₀, and constructing, based on the initial quasi-mean elements, an analytical perturbation model with a following mathematical expression:
  • assuming that α, i, Ω, ξ=e cos ω, η=−e sin ω, and λ=ω+M are first-type nonsingular osculating orbital elements of the object, where a represents the orbital semi-major axis, i represents the orbital inclination, Ω represents the orbital right ascension of ascending node, e represents the orbital eccentricity, ω represents the orbital argument of perigee, and M represents the orbital mean anomaly,

$\begin{array}{cc}a\left(t\right)={\stackrel{\_}{a}}_0+a_s^{\left(1\right)}\left(t\right)+a_s^{\left(2\right)}\left(t\right)& \left(1\right)\end{array}$ $\begin{array}{cc}i\left(t\right)={\stackrel{\_}{\iota }}_0+i_s^{\left(1\right)}\left(t\right)+i_s^{\left(2\right)}\left(t\right)& \left(2\right)\end{array}$ $\begin{array}{cc}\Omega \left(t\right)={\stackrel{\_}{\Omega }}_0+{\Omega }_1\left(t-T_0\right)+{\Omega }_s^{\left(1\right)}\left(t\right)+{\Omega }_s^{\left(2\right)}\left(t\right)& \left(3\right)\end{array}$ $\begin{array}{cc}\xi \left(t\right)={\stackrel{\_}{\xi }}_0\cos \left[{\omega }_1\left(t-T_0\right)\right]+{\stackrel{\_}{\eta }}_0\sin \left[{\omega }_1\left(t-T_0\right)\right]+{\xi }_s^{\left(1\right)}\left(t\right)+{\xi }_s^{\left(2\right)}\left(t\right)& \left(4\right)\end{array}$ $\begin{array}{cc}\eta \left(t\right)={\stackrel{\_}{\eta }}_0\cos \left[{\omega }_1\left(t-T_0\right)\right]-{\stackrel{\_}{\xi }}_0\sin \left[{\omega }_1\left(t-T_0\right)\right]+{\eta }_s^{\left(1\right)}\left(t\right)+{\eta }_s^{\left(2\right)}\left(t\right)& \left(5\right)\end{array}$ $\begin{array}{cc}\lambda \left(t\right)={\stackrel{\_}{\lambda }}_0+\left({\stackrel{\_}{n}}_0+{\lambda }_1\right)\left(t-T_0\right)+{\lambda }_s^{\left(1\right)}\left(t\right)+{\lambda }_s^{\left(2\right)}\left(t\right)& \left(6\right)\end{array}$

• • where in the above formulas, the left side represents osculating elements at a moment t, ā₀, l₀, Ω₀, ξ₀, η₀, and λ₀ represent initial quasi-mean elements at the moment T₀, n₀, =√{square root over (μ)}ā₀^−3/2 represents the angular velocity of mean motion of the object, μ represents the geocentric gravitational constant, Ω₁, ω₁, and λ₁ represent coefficients of first-order long-term variation items of corresponding elements, a_s⁽¹⁾, i_s⁽¹⁾, Ω_s⁽¹⁾, ξ_s⁽¹⁾, η_s⁽¹⁾, and λ_s⁽¹⁾ represent first-order short-period variation items of various elements, and a_s⁽²⁾, i_s⁽²⁾, Ω_s⁽²⁾, ξ_s⁽²⁾, η_s⁽²⁾, and λ_s⁽²⁾ represent second-order short-period variation items of the elements, which are reduced-order items related to rotation of the Earth.

Further, the forming an apparent path cluster by different apparent paths generated from the virtual objects during the current transit specifically includes:

• • taking the along-track error of the initial quasi-mean elements as a random variable uniformly distributed within [−τ, τ], discretizing the random variable to convert a probability problem into a certainty problem, and defining a following formula:

${\stackrel{\_}{\lambda }}_k={\stackrel{\_}{\lambda }}_0+\left({\stackrel{\_}{n}}_0+{\lambda }_1\right)k\mathrm{\Delta \tau }k=l,l-1,...,0,-1,...,-l+1,-l$

• • where Δτ represents a time increment, and assuming that ε>0 is a small value, a value of Δτ is as follows:

$\{\begin{array}{cc}\mathrm{\Delta \tau }=\epsilon & \tau \ge 2\epsilon \\ \mathrm{\Delta \tau }=\tau /2& \tau <2\epsilon \end{array}$

• • l is a positive integer, and is valued as follows:

$l=\mathrm{int}\left(\tau /\mathrm{\Delta \tau }\right)+1$

• • in the constructed analytical perturbation model, replacing λ₀ in the formula (6) with λ_k, and keeping other formulas unchanged to obtain a series of new perturbation models, which are denoted as a model _k, where each model _k corresponds to one virtual object, and a virtual object corresponding to ₀ is a theoretical object; and when k continuously changes from l to −l, generating, by each virtual object, a series of apparent paths to form an apparent path cluster {Γ_k|_l,−l}, where each apparent path Γ_k is uniquely determined by the model _k.

Further, the obtaining detection feature parameters of virtual objects corresponding to all valid apparent paths when the virtual objects rise to the search elevation specifically includes:

• • for one apparent path Γ_k in an apparent path cluster {Γ_k|_m,n}, where {Γ_k|_m,n} represents an apparent path cluster obtained by removing all the invisible apparent paths and further removing all the apparent paths with the unsatisfactory detectable arc length from {Γ_k|_l,−l}, performing search calculation forwards from a moment {tilde over (t)}_k of a closest approach point along a rising segment by using the dichotomous approach based on the model _k, generating, in each step, a position vector of the object relative to a True Equator and Mean Equinox of Epoch (TEMEE) coordinate system at a corresponding moment in the calculation process, calculating a corresponding detection elevation based on the obtained position vector and by considering the impact of the light travel time and an atmospheric refraction effect, determining a processing method for a next step based on a difference between the detection elevation and the search elevation, and finally obtaining detection feature parameters of the object when the object rises to the search elevation, including a detection moment, a detection range, and a detection azimuth; and repeating the above process for other apparent paths in the apparent path cluster to obtain detection feature parameters corresponding to all virtual objects.

Further, the generating a series of sets of search-specific guidance data specifically includes:

• • for one apparent path Γ_k in the apparent path cluster {Γ_k|_m,n}, where a detection moment, a detection range, and a detection azimuth of a corresponding virtual object when the virtual object rises to the search elevation are respectively t_k, ρ_k, and A′_k, and A′_k represents an azimuth after the continuation processing, taking t_k as a base point, and constructing, through cubic spline interpolation based on the known parameter ρ_k on the base point, an interpolation function in which a detection range ρ varies with a detection moment t; and taking A′_k as a base point, and constructing, through the cubic spline interpolation based on the known parameter t_k on the base point, an interpolation function in which the detection moment t varies with a detection azimuth A′, where the obtained interpolation functions are as follows:

$\begin{array}{cc}\rho =\rho \left(t\right)& \left(9\right)\end{array}$ $\begin{array}{cc}t=t\left(A^{\prime }\right)& \left(10\right)\end{array}$

• • a detection time scope (effective search period) when a real object rises to the search elevation is a closed interval [t_m, t_n]; a corresponding azimuth detection scope is a closed interval [A′_m, A′_n], which corresponds to an azimuth interval; and when the DLR system performs search along the search elevation, a change interval for a direction of the laser beam is [ψ_m, ψ_n], and a length of the interval is as follows:

$\mathrm{\Delta \psi }=❘{\psi }_m-{\psi }_n❘$

• • a length of the azimuth interval of the object on the search elevation is as follows:

$\Delta A=❘A_m^{\prime }-A_n^{\prime }❘$

• • obtaining, based on a differential geometry relationship, the length of the change interval for the direction of the laser beam:

$\mathrm{\Delta \psi }=\cosh \Delta A$

• • where h represents the search elevation of the DLR system;
  • evenly dividing the interval [ψ_m, ψ_n] into a plurality of subintervals based on an effective diameter w of the laser beam, where a quantity of the subintervals obtained through the division is as follows:

$N=\mathrm{int}\left[\frac{\mathrm{\Delta \psi }}{\left(1-\delta \right)w}\right]+1$

• • where δ represents a dimensionless scaling factor whose value is within an interval [0,1); and the azimuth interval and [ψ_m, ψ_n] are divided in a same manner, such that a length of an azimuth subinterval is as follows:

$\Delta A^{*}=\Delta A/N$

• • defining a symbol factor:

$\kappa =\{\begin{array}{cc}1& A_m^{\prime }<A_n^{\prime }\\ -1& A_m^{\prime }>A_n^{\prime }\end{array}$

• • where azimuth subintervals are sequentially [Ã*_v, Ã*_v+1], where v=1, 2, . . . , N, such that two endpoint values of each azimuth subinterval are obtained by recursion as follows:

${\tilde{A}}_1^{*}=A_m^{\prime }{\tilde{A}}_{v+1}^{*}={\tilde{A}}_v^{*}+\kappa \Delta A^{*}$

• • setting [t*_v, t*_v+1] to be a time subinterval corresponding to the azimuth subinterval [Ã*_v, Ã*_v+1], such that two endpoint values of the time subinterval are separately calculated according to the formula (10):

$t_v^{*}=t\left({\tilde{A}}_v^{*}\right)t_{v+1}^{*}=t\left({\tilde{A}}_{v+1}^{*}\right)$

• • calculating a center azimuth of the azimuth subinterval [Ã*_v, Ã*_v+1] according to a following formula:

${\tilde{A}}_{v+1/2}^{*}=\left({\tilde{A}}_v^{*}+{\tilde{A}}_{v+1}^{*}\right)/2$

• • calculating, according to the formula (10), a detection moment corresponding to the center azimuth:

$t_{v+1/2}^{*}=t\left({\tilde{A}}_{v+1/2}^{*}\right)$

• • calculating, according to the formula (9), a detection range corresponding to the center azimuth:

${\rho }_{v+1/2}^{*}=\rho \left({\tilde{A}}_{v+1/2}^{*}\right)$

• • restoring the center azimuth Ã*_v+1/2, which is a result obtained after the continuation processing, to an expression before the continuation processing, where a corresponding value A*_v+1/2 is as follows:

$A_{v+1/2}^{*}=\{\begin{array}{cc}{\tilde{A}}_{v+1/2}^{*}-2\pi & {\tilde{A}}_{v+1/2}^{*}\ge 2\pi \\ {\tilde{A}}_{v+1/2}^{*}+2\pi & {\tilde{A}}_{v+1/2}^{*}<0\\ {\tilde{A}}_{v+1/2}^{*}& 0\le {\tilde{A}}_{v+1/2}^{*}<2\pi \end{array}$

through the above processing process, generating a series of sets of search-specific guidance data A*_v+1/2, t*_v, t*_v+1, and ρ*_v+1/2, where v=1, 2, . . . , N, A*_v+1/2 represents an azimuth of the direction of the beam, t*_v and t*_v+1 respectively represent start and end dwell moments of the beam, and ρ*_v+1/2 represents a specified reference range of the object during the beam dwell.

Further, the performing, by the DLR system, multi-elevation search on the object based on the search-specific guidance data specifically includes:

• • generating a series of gradually increasing search elevations according to a principle that search elevations are generated in an ascending order and effective search periods of the search elevations do not intersect with each other, generating search-specific guidance data for each search elevation, and starting, by the DLR system based on the search-specific guidance data for the corresponding search elevation, search from a start search elevation and gradually switching to a highest search elevation, where
  • the series of gradually increasing search elevations are generated according to following steps:
  • 1) selecting a start search elevation h* between a minimum searchable elevation h_p and a maximum searchable elevation h_q; presetting a constant χ, where χ>0 represents adjustment time required for the DLR system to switch from one search elevation to another search elevation; and setting k=1, h_k=h*, and {tilde over (h)}=h*, where h_k represents a to-be-generated search elevation;
  • 2) calculating an effective search period [{tilde over (t)}_b, {tilde over (t)}_e] of the DLR system on the elevation {tilde over (h)}, and setting h_n={tilde over (h)} and h_x=h_q;
  • 3) setting h=(h_n+h_x)/2, calculating an effective search period [t_b, t_e] of the DLR system on the elevation h, and performing calculation according to ζ=t_b−{tilde over (t)}_e−χ;
  • 4) if |ζ|<ε_t, setting k=k+1, h_k=h, and {tilde over (h)}=h, and performing the step 2) to continue the calculation, where ε_t represents a given small value relative to χ; otherwise, performing next calculation;
  • 5) setting Δh=h_q−h; and if ζ<0 and Δ_h<ε_h, ending the calculation process, where En represents a given small value; otherwise, performing next calculation; and
  • 6) if ζ<0, setting h_n=h, and performing the step 3) to continue the calculation; otherwise, setting h_x=h, and performing the step 3) to continue the calculation; and
  • obtaining the series of gradually increasing search elevations h_k according to the above setting method.

Further, the procedure of calculating and correcting the along-track error of the initial quasi-mean elements based on a few sets of detection data obtained before switching to the tracking process, generating tracking-specific guidance data based on a set of corrected initial quasi-mean elements specifically includes:

• • setting one set of detection data obtained by the DLR system before switching to the tracking process as {circumflex over (t)} and {circumflex over (ρ)}, which is obtained during a certain dwell of the laser beam on a certain search elevation, where {circumflex over (t)} represents a detection moment, {circumflex over (ρ)} represents a corresponding detection range;
  • 1) setting an iterative initial value ϵ=0 of the along-track error, and using an atmospheric correction model to calculate a ranging delay Δρ caused by tropospheric refraction;
  • 2) performing calculation according to λ_ϵ=λ₀+(n₀+λ₁)ϵ, replacing λ₀ on the right side of the formula (6) with λ_ϵ in the analytical perturbation model ₀ to obtain a perturbation model _ϵ, and obtaining a calculated value ρ_c of the detection range corresponding to the detection moment {circumflex over (t)} and a time change rate {dot over (ρ)}_c of the calculated value based on the model _ϵ and in consideration of an impact of the light travel time; and
  • 3) correcting and updating the along-track error according to ϵ=ϵ+({circumflex over (ρ)}−Δρ−ρ_c)/{dot over (ρ)}_c, and then performing the step 2) to perform recalculation until along-track errors before and after the update meet certain convergence accuracy, to obtain one along-track error;
  • for each of L sets of detection data obtained before the DLR system switches to the tracking process, where L≥1, generating one along-track error ϵ according to the above calculation process, where a finally determined along-track error is as follows:

$\stackrel{\_}{ϵ}=\sum ϵ/L$

• • correcting the along-track error of the initial quasi-mean elements based on the along-track error ϵ, and generating the tracking-specific guidance data based on the corrected initial quasi-mean elements, which specifically includes:
  • performing calculation according to λ_ϵ=λ₀+(n₀+λ₁)ϵ, replacing λ₀ on the right side of the formula (6) with λ_ϵ in the analytical perturbation model ₀ to obtain a perturbation model _ϵ, and calculating, based on the model _ϵ, the impact of the light travel time, and an atmospheric refraction effect, detection azimuths, detection elevations, and detection ranges that correspond to a series of detection moments, where the detection azimuth and the detection elevation are used to determine a tracking direction of the DLR system, and the detection range is used to perform range-gate filtering on tracking echo signals of the DLR system.

The present disclosure further provides a computer-readable storage medium, where the computer-readable storage medium stores a computer program, and the computer program enables a computer to perform the above search and tracking method for full time-domain laser detection of space debris.

The present disclosure further provides an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor, where the processor, when executing the computer program, implements the above search and tracking method for full time-domain laser detection of space debris.

The present disclosure achieves following beneficial effects: The present disclosure adopts basic methods of celestial mechanics and astrometry, and constructs, based on precision orbital parameters of space debris, a new method for guiding a laser beam to search a debris object and switching to point tracking. Unlike various existing relevant laser detection methods of the space debris worldwide, the new method does not require initial guidance from an optical telescope with a large field of view when detecting the space debris, but can directly (blindly) detect and track the debris object using the laser beam, thereby breaking free from a time window limitation of an existing laser detection technology that relies on an optical device, and further achieving full time-domain laser detection of the space debris. Specific improvements of the present disclosure are as follows:

• • 1. An impact of the light travel time is considered during object search and tracking. Taking a low-orbit object with a range of 3000 kilometers as an example, a change to an apparent direction of the object caused by the light travel time is about 7 arc seconds, which is at a small magnitude relative to a beam diameter of a narrow-wave radar, and is at a same magnitude as a beam diameter of a DLR system. Therefore, during object search and tracking, the impact of the light travel time can be ignored for the narrow-wave radar, but cannot be ignored for the DLR system.
  • 2. A DLR beam is extremely narrow, and cumulative time for the object to enter the DLR beam in a constant elevation search process is very short, making it difficult to detect echo signals. The present disclosure proposes a time segment-based multi-elevation successive search scheme, which can effectively improve a detection success rate of the object.
  • 3. For further tracking after the object is detected, an existing patented technology (see patent document 1) generates an azimuth change rate and an elevation change rate that correspond to each dwelling beam in the search process. After detecting the object in the search process, the narrow-wave radar uses an azimuth change rate and an elevation change rate of a corresponding beam for linear propagation in real time to generate orbital predictions of the object. The orbital predictions generated in this manner is characterized by poor accuracy and fast attenuation. Although the orbital predictions can guide a radar beam for effective tracking for a short period, there is a significant technical risk involved in guiding an extremely narrow DLR beam. In this regard, the present disclosure proposes a novel prediction method. After detecting the object in the search process, the DLR system uses obtained detection data of the object in real time to generate an along-track error of a set of initial quasi-mean elements, corrects theoretical orbital elements based on the along-track error, and generates orbital predictions of the object based on the corrected orbital elements and an analytical perturbation model. Because an impact of an orbital dynamics factor is considered in the prediction, accuracy of the orbital predictions are significantly improved to guide the DLR beam to perform effective tracking for a long time. High-precision range information of the object in the orbital predictions can be used to significantly improve a detection probability of the DLR system for the echo signal, making DLR detection more effective. As effective tracking time increases, a requirement for real-time performance when the DLR system switches from a searching process to a tracking process is lowered, which is also conducive to detection and processing of the echo signal after the system detects the object.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 shows a geometric relationship indicating that a beam pointing angle differentially varies with an azimuth;

FIG. 2 is a processing flowchart of a search and tracking method for full time-domain laser detection of space debris;

FIG. 3 shows a setting of a multi-search elevation of AJISAI;

FIG. 4A shows that a pointing deviation (PD) and a range deviation (RD) of search of AJISAI vary with detection time, and FIG. 4B shows a locally enlarged view of a PD curve and its corresponding RD curve within a shaded block in FIG. 4A;

FIG. 5 shows that a PD and an RD of tracking of AJISAI vary with detection time;

FIG. 6 shows that a PD and an RD of prediction of AJISAI vary with detection time;

FIG. 7 shows a setting of a multi-search elevation of SWARM B;

FIG. 8A shows that a PD and an RD of search of SWARM B vary with detection time, and

FIG. 8B shows a locally enlarged view of a PD curve and its corresponding RD curve within a shaded block in FIG. 8A;

FIG. 9 shows that a PD and an RD of tracking of SWARM B vary with detection time; and

FIG. 10 shows that a PD and an RD of prediction of SWARM B vary with detection time.

DETAILED DESCRIPTION OF THE EMBODIMENTS

The technical solutions in the embodiments of the present disclosure are clearly and completely described below with reference to the accompanying drawings in the embodiments of the present disclosure. Apparently, the described embodiments are only some rather than all of the embodiments of the present disclosure. Based on the embodiments of the present disclosure, all other embodiments obtained by a person of ordinary skill in the art without creative efforts fall within the protection scope of the present disclosure.

Embodiment 1

As shown in FIG. 2, the present disclosure provides a search and tracking method for full time-domain laser detection of space debris, which is implemented during each transit of a debris object relative to a DLR system (only provided that constraints of a detection range and elevation are met). During implementation each time, it is required to obtain a set of latest precision orbital parameters determined based on data of precision radars and/or CCD telescopes of a space monitoring network, and start and end moments of a current transit of the object in advance. Search-specific guidance data of the DLR system is generated in advance based on the known information and in combination with estimation of a maximum along-track error of the orbital parameters of the object during the current transit. The DLR system performs constant elevation search on the object based on the search-specific guidance dat_q, obtains a plurality of pieces of detection data of the object after detecting the object during the search, determines an along-track error of the orbital parameters of the object based on the obtained detection dat_q, and corrects the orbital parameters of the object in real time based on the along-track error. Corrected orbital parameters are used to generate precise orbital predictions of the object to guide the DLR system to subsequently track and detect the object. The method specifically includes following steps.

Step 1: Modeling at an intermediate moment of the transit is performed.

Different types of orbital parameters can be converted to each other. It can be assumed that a set of known orbital parameters of the object are t_q, {right arrow over (r)}_q, {right arrow over ({dot over (r)})}_q, and t_q, where to represents an epoch moment of the orbital parameters, {right arrow over (r)}_q and {right arrow over (r)}_q respectively represent a position vector and a velocity vector of the object relative to an epoch geocentric inertial coordinate system, and ϵ_a represents an area-mass ratio of the object (which is a ballistic coefficient related to atmospheric drag perturbation). In addition, T_b and T_e are respectively set to start and end moments of a certain transit of the object. The intermediate moment of the current transit of the object is as follows:

$T_0=\frac{T_b+T_e}{2}$

The first-type nonsingular orbital elements are used as basic variables, and an analytical perturbation model with T₀ as an initial time is constructed. A model construction process is as follows:

• • 1) Perturbation propagation is performed by using a numerical method and a high-precision dynamic model based on the known orbital parameters t_q, {right arrow over (r)}_q, {right arrow over ({dot over (r)})}_q, and ta. Propagation is performed from the moment t_q to the moment T₀ to obtain position vector {right arrow over (r)}₀ and velocity vector {right arrow over ({dot over (r)})}₀ of the object relative to the epoch geocentric inertial coordinate system at the moment T₀. Conversion is performed on {right arrow over (r)}₀ and {right arrow over ({dot over (r)})}₀ to obtain a set of initial quasi-mean elements of the object at the moment T₀.
  • 2) Transit time of the object is relatively short, and there is a small perturbation change during the transit time. Therefore, the present disclosure constructs a simplified analytical perturbation model based on the initial quasi-mean elements at the moment T₀, which is used to describe orbital motion of the object during the current transit. The model only considers the first-order long-term items and the first-order short-period items of an orbital perturbation change, and a small quantity of second-order short-period items with a large amplitude due to rotation of the Earth.

A specific mathematical expression of the constructed model is as follows:

It is assumed that α, i, Ω, ξ=e cos ω, η=−e sin ω, and λ=ω+M are first-type nonsingular osculating orbital elements of the object, where a represents the orbital semi-major axis, i represents the orbital inclination, Ω represents the orbital right ascension of ascending node, e represents the orbital eccentricity, ω represents the orbital argument of perigee, and M represents the orbital mean anomaly. Therefore, following formulas are obtained:

$\begin{array}{cc}a\left(t\right)={\stackrel{\_}{a}}_0+a_s^{\left(1\right)}\left(t\right)+a_s^{\left(2\right)}\left(t\right)& \left(1\right)\end{array}$ $\begin{array}{cc}i\left(t\right)={\stackrel{\_}{\iota }}_0+i_s^{\left(1\right)}\left(t\right)+i_s^{\left(2\right)}\left(t\right)& \left(2\right)\end{array}$ $\begin{array}{cc}\Omega \left(t\right)={\stackrel{\_}{\Omega }}_0+{\Omega }_1\left(t-T_0\right)+{\Omega }_s^{\left(1\right)}\left(t\right)+{\Omega }_s^{\left(2\right)}\left(t\right)& \left(3\right)\end{array}$ $\begin{array}{cc}\xi \left(t\right)={\stackrel{\_}{\xi }}_0\cos \left[{\omega }_1\left(t-T_0\right)\right]+{\stackrel{\_}{\eta }}_0\sin \left[{\omega }_1\left(t-T_0\right)\right]+{\xi }_s^{\left(1\right)}\left(t\right)+{\xi }_s^{\left(2\right)}\left(t\right)& \left(4\right)\end{array}$ $\begin{array}{cc}\eta \left(t\right)={\stackrel{\_}{\eta }}_0\cos \left[{\omega }_1\left(t-T_0\right)\right]-{\stackrel{\_}{\xi }}_0\sin \left[{\omega }_1\left(t-T_0\right)\right]+{\eta }_s^{\left(1\right)}\left(t\right)+{\eta }_s^{\left(2\right)}\left(t\right)& \left(5\right)\end{array}$ $\begin{array}{cc}\lambda \left(t\right)={\stackrel{\_}{\lambda }}_0+\left({\stackrel{\_}{n}}_0+{\lambda }_1\right)\left(t-T_0\right)+{\lambda }_s^{\left(1\right)}\left(t\right)+{\lambda }_s^{\left(2\right)}\left(t\right)& \left(6\right)\end{array}$

In the above formulas, the left side represents osculating elements at moment t, ā₀, l₀, Ω₀, ξ₀, η₀, and λ₀ represent represent initial quasi-mean elements at the moment T₀, n₀, =√{square root over (μ)}ā₀^−3/2 represents an angular velocity of mean motion of the object, u represents the geocentric gravitational constant, Ω₁, ω₁, and λ₁ represent coefficients of first-order long-term perturbation items of corresponding elements, a_s⁽¹⁾, i_s⁽¹⁾, Ω_s⁽¹⁾, ξ_s⁽¹⁾, η_s⁽¹⁾, and λ_s⁽¹⁾ represent first-order short-period perturbation items of various elements, and a_s⁽²⁾, i_s⁽²⁾, Ω_s⁽²⁾, ξ_s⁽²⁾, η_s⁽²⁾, and λ_s⁽²⁾ represent second-order short-period perturbation items of the elements, which are reduced-order items related to the rotation of the Earth.

Step 2: A moment of a theoretical closest approach point is calculated.

A closest approach point of the object during the current transit corresponds to a certain point on the debris orbit. When the object reaches the point, the object has a maximum apparent elevation. The analytical perturbation model constructed in the step 1 determines a theoretical debris orbit, which corresponds to one theoretical closest approach point. The moment of the theoretical closest approach point is solved according to following two steps:

In a first step, moment t*₀ of an approximate closest approach point is solved according to a following formula:

$t_0^{*}=T_0+\frac{{\lambda }_s-{\lambda }_0}{\stackrel{.}{\lambda }-{\stackrel{.}{u}}_s}$

In the above formula, λ₀ represents a value of λ(t) at the moment T₀, which is the mean argument of latitude of the object at the moment T₀ and can be calculated according to the formula (6); {dot over (λ)}=n₀+λ₁ represents a long-term change rate of λ(t); u_s and λ_s respectively represent the true argument of latitude and the mean argument of latitude of a projection of an observation station on the orbit at the moment T₀; and {dot over (u)}_s represents a time change rate of u_s.

In a second step, the mean argument of latitude λ*₀ of the object at the moment t*₀ is first calculated according to the formula (6), and the true argument of latitude ud of the object at the moment t*₀ is calculated based on an elliptical motion relationship and λ*₀. Then, u*₀ is taken as an initial value to perform iterative solving to obtain the true argument of latitude to ũ₀ corresponding to moment {tilde over (t)}₀ of an accurate closest approach point. An equation for the iterative solving is as follows:

$\sin \left({\tilde{u}}_0-u_0^{*}\right)=\frac{\left({\tilde{\xi }}_0\sin {\tilde{u}}_0+{\tilde{\eta }}_0\cos {\tilde{u}}_0\right){\sin }^2\theta }{\left(1-{\tilde{\eta }}_0\sin {\tilde{u}}_0+{\tilde{\xi }}_0\cos {\tilde{u}}_0\right)\cos {\theta }_0\left(r/R-\cos \theta \right)}$

In the above equation, {tilde over (ξ)}₀ and {tilde over (η)}₀ respectively represent values of ξ(t) and η(t) at the moment {tilde over (t)}₀, which can be respectively calculated according to the formulas (4) and (5); θ represents a geocentric angle of the object and the observation station at the moment {tilde over (t)}₀; θ₀ represents a latitude of the observation station with reference to orbit plane at the moment t*₀; r represents a geocentric distance of the object at the moment {tilde over (t)}₀; and R represents a geocentric distance of the observation station. Finally, the mean argument of latitude {tilde over (λ)}₀ of the object at the moment {tilde over (t)}₀ is calculated based on the elliptical motion relationship and ũ₀and the moment of the accurate closest approach point is calculated according to a following formula:

${\tilde{t}}_0=t_0^{*}+\frac{{\tilde{\lambda }}_0-{\lambda }_0^{*}}{\stackrel{.}{\lambda }}$

Step 3: The along-track error is discretized and an apparent path cluster is generated.

Based on accuracy of the known orbital parameters t_q, {right arrow over (r)}_q, {right arrow over ({dot over (r)})}_q, ϵ_a, and prediction duration, a maximum along-track error of the initial quasi-mean elements at the moment T₀ is estimated. A corresponding estimated value is a time quantum whose absolute value is set to τ. The along-track error of the initial quasi-mean elements is considered as a random variable uniformly distributed within closed interval [−τ, τ], and the random variable is discretized to convert a probability problem into a certainty problem. Therefore, a following formula is defined:

${\stackrel{\_}{\lambda }}_k={\stackrel{\_}{\lambda }}_0+\left({\stackrel{\_}{n}}_0+{\lambda }_1\right)k\mathrm{\Delta \tau }k=l,l-1,...,0,-1,...,-l+1,-l$

In the above formula, Δτ represents a time increment, and assuming that ε>0 is a small value, a value of Δτ is as follows:

$\{\begin{array}{cc}\mathrm{\Delta \tau }=\epsilon & \tau \ge 2\epsilon \\ \mathrm{\Delta \tau }=\tau /2& \tau <2\epsilon \end{array}$

• • l is a positive integer, and is valued as follows:

$l=\mathrm{int}\left(\tau /\mathrm{\Delta \tau }\right)+1$

In the constructed analytical perturbation model, λ₀ on the right side of the formula (6) is replaced with λ_k, and other formulas keep unchanged to obtain a new perturbation model, which is denoted as model _k. Each model _k determines one virtual object, and a virtual object corresponding to ₀ is a theoretical object. Each virtual object moves along the same theoretical orbit, with only different along-track positions. Due to the rotation of the Earth, from a ground station, virtual objects at different positions generate different apparent motion tracks during the current transit of the object. When k continuously changes from l to −l, the corresponding virtual object generates a series of apparent paths to form an apparent path cluster {Γ_k|_l,−l}, where each apparent path Γ_k is uniquely determined by the model _k.

Step 4: Detectability of the closest approach point is constrained.

The closest approach point is a point on the apparent path with the best detection condition for the DLR system. If the closest approach point is not detectable, the entire apparent path cannot be detected, and such an apparent path cannot be considered. Therefore, it is necessary to confirm detectability of each apparent path in the apparent path cluster {Γ_k|_l,−l} and remove all undetectable apparent paths.

Firstly, a moment of a closest approach point of each apparent path is calculated as follows:

• • 1) Theoretical apparent path Γ₀ is considered, and its detectability has been confirmed in prediction calculation on the transit. A moment of a corresponding closest approach point has been calculated in the step 2.
  • 2) Two adjacent apparent paths Γ_k and Γ_k+1 are considered, a moment of an accurate closest approach point of Γ_k is used as a moment of an approximate closest approach point of Γ_k+1, and a second step of the step 2 is directly preformed to enter a solving process. In the solving process, model _k+1 is used for relevant calculation to obtain a moment of an accurate closest approach point of Γ_k+1. k is considered to continuously change from 0 to l−1. The above processing method is used to gradually obtain moments of closest approach points of apparent paths Γ₁, Γ₂, and . . . , Γ_l.
  • 3) Two adjacent apparent paths Γ_k and Γ_k−1 are considered, the moment of the accurate closest approach point of Γ_k is used as a moment of an approximate closest approach point of Γ_k−1, and the second step of the step 2 is directly preformed to enter the solving process. In the solving process, model _k−1 is used for relevant calculation to obtain a moment of an accurate closest approach point of Γ_k−1, k is considered to continuously change from 0 to −l +1. The above processing method is used to gradually obtain moments of closest approach points of apparent paths Γ₋₁, Γ₋₂, and . . . , Γ_−l.

So far, moments of closest approach points of all apparent paths are calculated.

For one apparent path Γ_k in the apparent path cluster {Γ_k|_l,−l}, moment {tilde over (t)}_k of a closest approach point of the one apparent path is substituted into the model _k to obtain osculating elements at the moment of the closest approach point. A range and an elevation of the closest approach point are calculated based on the osculating elements, and are compared with detection constraints of the DLR system to confirm detectability of one apparent path. The detection constraints of the DLR system include thresholds of the detection range and elevation. After detectability of each apparent path in the apparent path cluster {Γ_k|_l,−l} is confirmed, all the undetectable apparent paths are removed, and the remaining apparent paths are still consecutive, forming apparent path cluster {Γ_k|_i,j}, where i≥0≥j.

Step 5: Detection effectiveness is constrained.

A detectable arc length of the DLR system is a necessary condition for the DLR system to obtain high-quality detection data. Therefore, it is necessary to further screen each apparent path in the apparent path cluster {Γ_k|_i,j} and remove all apparent paths with unsatisfactory detectable arc lengths to effectively obtain DLR data and improve detection efficiency of the system.

Taking the moment of the closest approach point as a boundary, apparent motion of the object is a process of first increasing the elevation and then decreasing the elevation. An apparent path segment before the moment of the closest approach point is an ascending segment, and an apparent path segment after the moment of the closest approach point is a descending segment. In order to ensure that the DLR system can still have a long detectable arc segment after searching and discovering the object, the present disclosure is intended to search an object in the ascending segment, and require that detectable duration of an object in the descending segment should not be less than a given threshold. A specific processing process is as follows:

For one apparent path Γ_k in the apparent path cluster {Γ_k|_i,j}, based on the model _k and in combination with the detection constraints of the DLR system, search calculation is carried out backwards along the descending segment from the moment {tilde over (t)}_k of the closest approach point at an appropriate step until a moment specified by a threshold. If an undetectable target point is found during the search, the one apparent path is removed from the apparent path cluster. The above processing process is repeated for other apparent paths in the apparent path cluster, such that all the apparent paths with the unsatisfactory detectable arc lengths are removed. If there is no remaining apparent path, it indicates that the DLR system cannot obtain valid detection data during the current transit of the object, and subsequent processing is terminated. Otherwise, the remaining apparent paths are still consecutive, forming apparent path cluster {Γ_k|_m,n}, where m≥n.

Step 6: A range of a search elevation is determined.

For one apparent path Γ_k in the apparent path cluster {Γ_k|_m,n}, a maximum detectable elevation of the one apparent path is an elevation of a corresponding closest approach point, which has been given in the step 4. A minimum detectable elevation of the one apparent path needs to be calculated. Search calculation can be performed forwards along the ascending segment from the moment {tilde over (t)}_k of the closest approach point at an appropriate initial step by using a dichotomous approach based on the model _k. During the calculation, in each step, detectability can be determined based on the detection constraints of the DLR system to determine a processing method for a next step. Finally, a minimum detectable elevation meeting a certain accuracy requirement is given. The above processing process is repeated for other apparent paths in the apparent path cluster, such that maximum detectable elevations and minimum detectable elevations of all apparent paths are obtained.

A minimum value is taken from the maximum detectable elevations of all the apparent paths, and set to h_q. A maximum value is taken from the minimum detectable elevations of all the apparent paths, and set to h_p. Then an elevation scope suitable for the DLR system to perform constant elevation search is determined based on closed interval [h_p, h_q], where h_q≥h_p. An available search elevation (ASE) of the DLR system can be obtained according to a following formula:

$h=h_p+\beta \left(h_q-h_p\right)$

In the above formula, β represents a scaling factor for taking a value in closed interval [0, 1].

Step 7: Feature parameters on the search elevation are calculated.

A space object travels to a target point along its orbit at moment t′. A position vector of the target point relative to the geocentric inertial coordinate system is {right arrow over (r)}. The laser beam emitted by the DLR system at the moment t hits the target point, where t represents the detection moment of the DLR system. The corresponding detection distance, azimuth, and elevation of the detection moment can be calculated based on position vector {right arrow over (r)}_b of the target point relative to the Earth-fixed coordinate system and coordinates of the ground station. Ty can be obtained through following coordinate conversion:

$\begin{array}{cc}{\overrightarrow{r}}_b=G\left(t\right)\overrightarrow{r}& \left(7\right)\end{array}$

As described above, G represents a conversion matrix from the geocentric inertial coordinate system to the Earth-fixed coordinate system, which is a time function. Due to a limited speed of light, a light wave requires a propagation time from the ground station to the target point, which is referred to as the light travel time. If an impact of the light travel time is ignored, it can be set that t=t′. Then {right arrow over (r)}_b can be directly obtained according to the formula (7), and then a detection range, azimuth, and elevation of the target point can be given. If the impact of the light travel time is considered, a following formula is obtained:

$\begin{array}{cc}t=t^{\prime }-\rho /c& \left(8\right)\end{array}$

In the above formula, ρ represents the detection range of the target point, c represents the speed of light, and t and {right arrow over (r)}_b are interdependent, which are not known and can only be iteratively processed. During the iteration, initial condition t=t′ is set, {right arrow over (r)}_b is calculated according to the formula (7), and a new value of t is obtained according to the formula (8) until t meets specified convergence accuracy. Finally, an accurate detection moment and its corresponding detection range, azimuth, and elevation are given.

For one apparent path Γ_k in the apparent path cluster {Γ_k|_m,n}, search calculation is performed forwards along the ascending segment from the moment t_k of the closest approach point at an appropriate initial step by using the dichotomous approach based on the model _k. During the calculation, in each step, a position vector of the object relative to a TEMEE coordinate system at the corresponding moment is generated first. Then, a corresponding detection elevation is calculated based on the obtained position vector and by considering the impact of the light travel time and an atmospheric refraction effect, a processing method for a next step is determined based on a difference between the detection elevation and the search elevation, and finally detection feature parameters meeting a certain accuracy requirement are given, including a detection moment, range, and azimuth when the object rises to the search elevation. The above processing process is repeated for other apparent paths in the apparent path cluster, such that detection feature parameters corresponding to all virtual objects are obtained.

Step 8: Continuation processing is performed on the azimuth.

For one apparent path Γ_k in the apparent path cluster {Γ_k|_m,n}, an azimuth of a corresponding virtual object when the virtual object rises to the search elevation is A_k, where a value of A_k is within interval [0,2π). The value of A_k is redefined, and a defined value of A_k is A′_k.

k is considered to continuously change from m to n, and A′_m=A_m and ΔA_k=A_k−A_k+1 are set, where m−1≥k≥n. Subsequent A′_k is calculated based on following three independent cases:

• • 1) If there is integer p that meets m−1≥p≥n, and ΔA_p<−π, a following formula is obtained:

$A_k^{\prime }=\{\begin{array}{cc}{A}_k& k>p\\ A_k+2\pi & k\le p\end{array}$

• • 2) If there is integer p that meets m−1≥p≥n, and ΔA_p>π, a following formula is obtained:

$A_k^{\prime }=\{\begin{array}{cc}{A}_k& k>p\\ A_k-2\pi & k\le p\end{array}$

• • 3) If there is no integer p meeting the above two cases, a following formula is obtained:

$\begin{array}{cc}{A}_k^{\prime }=A_k& m-1\ge k\ge n\end{array}$

A value of each A′_k obtained above is no longer limited to the interval [0,2λ), but is consecutive while varying with k. A′_k and A_k correspond to a same physical orientation, and A′_k can be used to replace A_k as the azimuth.

Step 9: The search-specific guidance data is generated.

For one apparent path Γ_k in the apparent path cluster {Γ_k|_m,n}, a detection moment, range, and azimuth when a corresponding virtual object rises to the search elevation are respectively set to t_k, ρ_k, and A′_k, where A′_k represents an azimuth after the continuation processing. k is considered to continuously change from m to n, and t_k strictly monotonically increases with k. Therefore, t_k can be used as a base point. An interpolation function in which the detection range ρ varies with the detection moment t is constructed through cubic spline interpolation based on the known parameter ρ_k on the base point. A′_k strictly monotonically increases or decreases with k. Therefore, A′_k can be used as a base point, and an interpolation function in which the detection moment t varies with the detection azimuth A′ is constructed through the cubic spline interpolation based on the known parameter t_k on the base point. The obtained interpolation functions are as follows:

$\begin{array}{cc}\rho =\rho \left(t\right)& \left(9\right)\end{array}$ $\begin{array}{cc}t=t\left(A^{\prime }\right)& \left(10\right)\end{array}$

A detection time scope when a real object rises to the search elevation is closed interval [t_m, t_n], and a corresponding azimuth detection scope is closed interval [A′_m, A′_n], which corresponds to an azimuth interval. When the DLR system performs search along the search elevation, there is also a change interval for a direction of the laser beam, which is denoted as [ψ_m, ψ_n], and a length of the interval is as follows:

$\mathrm{\Delta \psi }=❘{\psi }_m-{\psi }_n❘$

A length of the azimuth interval of the object on the search elevation is as follows:

$\Delta A=❘A_m^{\prime }-A_n^{\prime }❘$

As shown in FIG. 1, a following differential geometric relationship is used:

$d\psi =\cosh \mathrm{dA}$ ${\int }_{{\psi }_m}^{{\psi }_n}d\psi ={\int }_{A_m^{\prime }}^{A_n^{\prime }}\cosh \mathrm{dA}$

The length of the change interval for the direction of the laser beam can be obtained:

$\mathrm{\Delta \psi }=\cosh \Delta A$

An effective diameter of the laser beam is set to w. The interval [ψ_m, ψ_n] is evenly divided into a plurality of subintervals. A quantity of the subintervals obtained through the division is as follows:

$N=\mathrm{int}\left[\frac{\mathrm{\Delta \psi }}{\left(1-\delta \right)w}\right]+1$

In the above formula, δ represents a dimensionless scaling factor that takes a value within interval [0,1), and a length of the subinterval [ψ_m, ψ_n] obtained through the division based on N does not exceed the effective diameter of the laser beam. The azimuth interval and [ψ_m, ψ_n] are divided in a same manner, such that a length of an azimuth subinterval is as follows:

$\Delta A^{*}=\Delta A/N$

A symbol factor is defined:

$\kappa =\{\begin{array}{cc}1& A_m^{\prime }<A_n^{\prime }\\ -1& A_m^{\prime }>A_n^{\prime }\end{array}$

Assuming that azimuth subintervals are sequentially [Ã*_v, Ã*_v+1], where v=1, 2, . . . , N, two endpoint values of each azimuth subinterval are obtained by recursion as follows:

${\tilde{A}}_1^{*}=A_m^{\prime }$ ${\tilde{A}}_{v+1}^{*}={\tilde{A}}_v^{*}+\kappa \Delta A^{*}$

• • [t*_v, t*_v+1] is set to be a time subinterval corresponding to the azimuth subinterval [Ã*_v, Ã*_v+1], such that two endpoint values of the time subinterval are separately calculated according to the formula (10):

$t_v^{*}=t\left({\tilde{A}}_v^{*}\right)$ $t_{v+1}^{*}=t\left({\tilde{A}}_{v+1}^{*}\right)$

A center azimuth of the azimuth subinterval [Ã*_v, Ã*_v+1] is calculated according to a following formula:

${\tilde{A}}_{v+1/2}^{*}=\left({\tilde{A}}_v^{*}+{\tilde{A}}_{v+1}^{*}\right)/2$

A detection moment corresponding to the center azimuth is calculated according to the formula (10):

${\tau }_{v+1/2}^{*}=\tau \left({\tilde{A}}_{v+1/2}^{*}\right)$

A detection range corresponding to the center azimuth is calculated according to the formula (9):

${\rho }_{v+1/2}^{*}=\rho \left(t_{v+1/2}^{*}\right)$

The center azimuth A*_v+1/2 is actually a result obtained after the continuation processing, and is restored to a common expressive method. Assuming that a corresponding value is Ã*_v+1/2, a following formula is obtained:

$A_{v+1/2}^{*}=\{\begin{array}{cc}{\tilde{A}}_{v+1/2}^{*}-2\pi & {\tilde{A}}_{v+1/2}^{*}\ge 2\pi \\ {\tilde{A}}_{v+1/2}^{*}+2\pi & {\tilde{A}}_{v+1/2}^{*}<0\\ {\tilde{A}}_{v+1/2}^{*}& 0\le {\tilde{A}}_{v+1/2}^{*}<2\pi \end{array}$

Through the above processing process, a series of sets of search-specific guidance data A*_v+1/2, t*_v, t*_v+1, and ρ*_v+1/2 is generated, where v=1,2, ... , N. Each set of guidance data corresponds to one dwell of the laser beam on the search elevation, where A*_v+1/2 represents an azimuth of the direction of the beam, t*_v and t*_v+1 respectively represent start and end dwell moments of the beam, and ρ*_v+1/2 represents a specified reference range of the object during the beam dwell, which is used to perform range-gate filtering on echo signals. When each set of guidance data changes from v=1 to v=N, the corresponding beam is considered to dwell in different orientations on the search elevation in chronological order, forming a method and process of the constant elevation search by the DLR system. A start time of a first dwell of the beam is t*₁=t_m. For each of other beams, a start time of a current dwell of the beam is an end time of a previous dwell of the beam. An end time of a last dwell of the beam is t*_N+1=t_n, and [t_m, t_n] is an effective search period of the DLR system on the search elevation.

The above is only a general processing process, and there is a special case that needs to be explained to ensure preciseness and completeness of the technical solution:

After the processing in the step 5, if there is only one apparent path in the apparent path cluster {Γ_k|_m,n}, a detection moment when a virtual object on the one apparent path rises to the search elevation is taken as a center moment, and the center moment is separately subtracted and added by Δτ/2 as the start and end dwell moments of the beam. There is a total of one dwell for the beam, and the azimuth of the direction of the beam and the reference range set for the object during the beam dwell are taken as corresponding detection feature parameters of the virtual object on the one apparent path.

Step 10: A plurality of search elevations are set.

Multi-elevation search is an effective means to improve an object detection success rate of the DLR system. Based on a motion characteristic of the object in the ascending segment, the DLR system should carry out the multi-elevation search gradually from a low elevation to a high elevation. In addition, effective search periods for the plurality of elevations should not intersect with each other to avoid a time conflict between different search processes. Following the above consideration, a method for setting the plurality of search elevations is given as follows:

• • 1) Start search elevation h* between minimum searchable elevation h_p and maximum searchable elevation h_q is selected based on a technical condition of the DLR system, and generally, h*≥15°. In addition, constant χ is preset, where χ>0 represents adjustment time required for the DLR system to switch from one search elevation to another search elevation. It is set that k=1, h_k=h*, and {tilde over (h)}=h*.
  • 2) Effective search period [{tilde over (t)}_b, {tilde over (t)}_e] of the DLR system on the elevation {tilde over (h)} is obtained through relevant calculation in the step 7, and h_n={tilde over (h)} and h_x=h are set.
  • 3) It is set that h=(h_n+h_x)/2, effective search period [t_b, t_e] of the DLR system on the elevation h is calculated through relevant calculation in the step 7, and calculation is performed according to ζ=t_b−{tilde over (t)}_e−χ.
  • 4) If |ζ|<ε_t, k=k+1, h_k=h are set, and the step 2) is performed to continue the calculation, where ε_t represents a given small value relative to x; otherwise, next calculation is performed.
  • 5) It is set that Δh=h_q−h; and if ζ<0 and Δ_h<ε_h, the calculation process ends, where En represents a given small value; otherwise, next calculation is performed.
  • 6) If ζ<0, h_n=h is set, and the step 3) is performed to continue the calculation; otherwise, h_x=h is set, and the step 3) is performed to continue the calculation.

By using the above setting method, a series of gradually increasing search elevations h_k are obtained, where k≥1, and a maximum quantity of effective search elevations are set. Search-specific guidance data for each elevation can be generated by repeatedly performing the steps 7, 8, and 9. Based on search-specific guidance data for a corresponding elevation, the DLR system starts search on a start elevation and can gradually switch to perform search on a higher elevation, until search on a highest elevation is completed.

Step 11: Tracking-specific guidance data is generated.

After discovering the object for the first time during the search, the DLR system can switch to a tracking process at any time. It is assumed that the DLR system has obtained a few sets of detection data before switching to the tracking process, and one set of detection data is {circumflex over (t)} and {circumflex over (ρ)}, which is obtained during a certain dwell of the laser beam on a certain search elevation, where {circumflex over (t)} represents the detection moment, and {circumflex over (ρ)} represents the corresponding detection range. In order to generate high-precision guidance data suitable for DLR tracking, it is necessary to correct the along-track error for the initial quasi-mean elements. One along-track error ϵ can be obtained through iterative calculation based on the one set of detection dat_q, and a calculation process is as follows:

• • 1) Iterative initial value ϵ=0 is set, and an appropriate atmospheric correction model is used to calculate ranging delay Δρ caused by tropospheric refraction.
  • 2) Calculation is performed according to λ_ϵ=λ₀+(n₀+λ₁)ϵ, λ₀ on the right side of the formula (6) is replaced with λ_ϵ and other formulas keep unchanged in the analytical perturbation model ₀ to obtain perturbation model _ϵ, and calculated value ρ_c of the detection range corresponding to the detection moment {circumflex over (t)} and time change rate {dot over (ρ)}_c of the calculated value are obtained based on the model_ϵ and in consideration of an impact of the light travel time.
  • 3) The along-track error is corrected and updated according to ϵ=ϵ+({circumflex over (ρ)}−Δρ−ρ_c)/{dot over (ρ)}_c, and then the step 2) is performed for recalculation until along-track errors before and after the update meet certain convergence accuracy, to obtain one along-track error.

Based on each of L sets of detection data obtained before the DLR system switches to the tracking process, where L≥1, one along-track error e is generated according to the above calculation process. A finally determined along-track error is as follows:

$\stackrel{\_}{ϵ}=\sum ϵ/L$

The along-track error is corrected for the initial quasi-mean elements based on error ϵ, and the tracking-specific guidance data is generated based on a set of corrected initial quasi-mean elements. A specific method includes: performing calculation according to λ_ϵ=λ₀+(n₀+λ₁)ϵ, replacing λ₀ on the right side of the formula (6) with λ_ϵand keeping other formulas unchanged in the analytical perturbation model ₀ to obtain perturbation model _ϵ, and calculating, based on the model _ϵ, the impact of the light travel time, and the atmospheric refraction effect, detection azimuths, detection elevations, and detection ranges that correspond to a series of detection moments. The detection azimuth and the detection elevation are used to determine a tracking direction of the DLR system, and the detection range is used to perform range-gate filtering on tracking echo signals of the DLR system.

The present disclosure always implements the processing for the along-track error. The estimated maximum along-track error is discretized to form a constant elevation search process of the DLR system. The DLR system searches and discovers the object during the search, and the detection data obtained by the system is used as posterior information, to determine the along-track error in real time. The along-track error is corrected for theoretical orbital elements to generate high-precision tracking-specific guidance data.

Embodiment 2

This embodiment selects AJISAI (NASA catalog number 16908) as an experimental object, which has an orbital height of 1,400 kilometers and is a laser calibration satellite. A high-precision CPF ephemeris of the object can be used as an experimental comparison standard. The experiment selects a hypothetical DLR site in China and assumes that an effective beam diameter of a DLR system is 20 arc seconds. When the object enters an effective DLR beam, an included angle between a beam direction and a direction of the object should be less than 10 arc seconds. In the experiment, one set of orbital parameters of the object relative to the J2000 geocentric inertial coordinate system is generated based on detection data obtained by a precision tracking radar in China, namely t_q, {right arrow over (r)}_q, _q, and ϵ_a, as shown in Table 1:

TABLE 1
One set of orbital parameters of the AJISAI relative
to the J2000 geocentric inertial coordinate system
tq(Beijing time)	Feb. 25, 2021 07:50:56.091
{right arrow over (r)}q	xq	1059237.065	q	{dot over (x)}q	−6942.578
(m)	yq	4948397.066	(m/s)	{dot over (y)}q	1557.054
	zq	6023858.754		żq	−55.253
ϵa(m2/kg)	0.00010

The experiment is based on the one set of orbital parameters and specified detection constraints of the DLR system (a maximum detection range is 5,000 kilometers, and a minimum detection elevation is 5°) , and a certain transit of the object is determined through high-precision perturbation calculation. The transit is a transit after 19 days starting from the epoch moment of the orbital parameters of the object. During the transit, search, and tracking experiments of the object are carried out.

1) Intermediate moment T₀ of the current transit is calculated first. In addition, perturbation propagation is performed by using a numerical method and a high-precision dynamic model based on the known orbital parameters t_q, {right arrow over (r)}_q, {right arrow over ({dot over (r)})}_q, and ϵ_a. Propagation is performed from the moment t_q to the moment T₀ to obtain position vector {right arrow over (r)}₀ and velocity vector {right arrow over ({dot over (r)})}₀ of the object relative to the J2000 geocentric inertial coordinate system at the moment T₀. Further, {right arrow over (r)}₀ and {right arrow over ({dot over (r)})}₀ are converted into a set of initial quasi-mean elements of the object at the moment T₀, as shown in Table 2.

TABLE 2
Initial quasi-mean elements of AJISAI at the intermediate
moment during the transit 10:15:25.450
T0 (Beijing time)	Mar. 16, 2021 10:15:25.450
ā0(m)	7866388.521	ι 0	50°.0079
Ω 0	288°.1634	ξ 0	0.0010384276
η 0	−0.0002674244	λ 0	32°.2032

Based on the initial quasi-mean elements, analytical perturbation model ₀ describing theoretical orbital motion of the object during the current transit is determined according to formulas (1) to (6).

2) Firstly, a geocentric distance and a geocentric longitude and latitude of an observation station are calculated based on geodetic coordinates of the DLR system. Then, based on the analytical perturbation model ₀, two-step calculation is performed. In a first step, moment t*₀ of an approximate closest approach point is calculated based on T₀. In a second step, moment {tilde over (t)}₀ of an accurate closest approach point is obtained through iterative calculation based on t*₀. The moment of the accurate closest approach point is a moment of a theoretical closest approach point of the object during the current transit.

3) It is estimated that absolute value τ of a maximum along-track error of the initial quasi-mean elements is 3 seconds, and ε is set to 0.5 seconds. Therefore, Δτ is 0.5 seconds, and it is obtained through calculation that l=7. Then, calculation is performed according to a following formula:

${\overline{\lambda }}_k={\overline{\lambda }}_0+\left({\overline{n}}_0+{\lambda }_1\right)k\Delta \tau$ $k=7,6,\dots ,0,\dots ,-6,-7$

In the analytical perturbation model ₀, λ_k is used to replace λ₀ in turn to obtain a series of perturbation models _k. Each model _k describes motion of a virtual object, and apparent path cluster {Γ_k|_{7, −7}} is formed by apparent motion tracks of all virtual objects during the current transit.

4) Firstly, a moment of a closest approach point of each apparent path in the apparent path cluster {Γ_k|_{7, −7}} is calculated. Then, a range and an elevation of the closest approach point of each apparent path are obtained. The visibility of the closest approach point of each apparent path is confirmed, and no undetectable apparent path is found. All undetectable apparent paths are removed, and an apparent path cluster formed by the remaining apparent paths is still {Γ_k|_{7, −7}}.

5) Detectable duration of an object in a descending segment is limited to not less than 10 seconds. For one apparent path Γ_k in the apparent path cluster {Γ_k|_{7, −7}}, based on the model _k, search calculation is carried out backwards along the descending segment from moment {tilde over (t)}_k of a closest approach point of the one apparent path at a step of 10 seconds. In the calculation, only one step is performed to obtain a target range and elevation on a corresponding calculation point, then are compared with the detection constraints of the DLR system. If not matched, the one apparent path is removed from the apparent path cluster; otherwise, the one apparent path is retained. For other apparent paths in the apparent path cluster {Γ_k|_{7, −7}}, the above processing process is repeated. An apparent path cluster formed by the finally remaining apparent paths is still {Γ_k|_{7, −7}}.

6) For one apparent path Γ_k in the apparent path cluster {Γ_k|_{7, −7}}, a maximum detectable elevation of the one apparent path is an elevation of a closest approach point of the one apparent path, and a minimum detectable elevation of the one apparent path is obtained through search in an ascending segment of the one apparent path by using a dichotomous approach. During search calculation, an initial step of the dichotomous approach is set to 30 seconds, and convergence accuracy is set to 10⁻³ seconds. For other apparent paths in the apparent path cluster {Γ_k|_{7, −7}}, the above processing process is repeated to finally obtain maximum detectable elevations and minimum detectable elevations of all apparent paths. Maximum value h_p is taken from the minimum detectable elevations of all the apparent paths, and minimum value h_q is taken from the maximum detectable elevations of all the apparent paths, such that following values are obtained:

$h_p=5°·5000$ $h_q=79°·9082$

Closed interval [h_p, h_q] is an elevation scope suitable for the DLR system to perform constant elevation search.

7) Starting search elevation h*=15° is selected from [h_p, h_q]. It is set that χ=10 seconds, ε_t=0.1 seconds, and ε_h=2°. A relevant method in the step 7 of Embodiment 1 is used to calculate effective search period [t_b, t_e] for elevation {tilde over (h)} and effective search period [t_b, t_e] for elevation h. Finally, 26 search elevations with non-intersecting effective search periods are obtained. Based on one obtained search elevation, the steps 7, 8, and 9 of Embodiment 1 are performed in sequence to generate search-specific guidance data for the search elevation. Based on other obtained search elevations, the above processing process is repeated to ultimately generate search-specific guidance data for all the 26 search elevations. In the above processing process, when the step 7 of Embodiment 1 is performed for calculation, the initial step of the dichotomous approach is set to 30 seconds, and the convergence accuracy is set to 10⁻⁴ degrees. An atmospheric refraction effect contained in detection feature parameters is provided by the Hopfield model based on a ratio of 90%. When the step 9 of Embodiment 1 is performed for calculation, scaling factor δ is set to 0.1.

FIG. 3 shows search elevations whose degrees are in an ascending order and their corresponding effective search periods in chronological order from left to right. As shown in FIG. 3, one of the search elevations is 45.5797 degrees, and a corresponding effective search period starts from 10:12:6.893 on Mar. 16, 2021 and ends at 10:12:14.164 on Mar. 16, 2021. A total of 19 sets of search-specific guidance data are generated for the elevation. The DLR system chronologically enables, based on each set of search-specific guidance dat_q, the beam to dwell, to form a constant elevation search process.

For any detection moment during the current transit of the object, a target point at the moment can be generated by combining CPF ephemeris interpolation and the light travel time calculation. The atmospheric refraction effect when the DLR system detects the target point is considered, and the Hopfield model is used to calculate a refractive value to obtain one detection point corresponding to the target point, which is referred to as a standard point of the object at the moment. For any detection moment in the above constant elevation search process, a PD of the search of the object at the moment is defined as an included angle between a direction of a dwelling beam of the DLR system at the moment and a direction of a standard point of the object at the same moment, and an RD of the search of the object at the moment is defined as an absolute value of a difference between a reference range that is of the object and specified for the dwelling beam of the DLR system at the moment and a slant range of the standard point of the object at the same moment. A situation of detecting the object by the DLR system in the above constant elevation search process can be determined based on PDs of the search of the object at different detection moments. FIG. 4A and FIG. 4B are coordinate graphs with two Y-axes, which respectively show changes to the PD of search of the object and the RD of the search of the object over time in the above constant elevation search process. FIG. 4A shows change curves of the PD of search and the RD of the search over time in the entire search process. From FIG. 4A, it can be seen that the object enters the dwelling beam of the DLR system during a very short period in the search process, and the RD of the search during this period is also very small, indicating that the DLR system can detect the object in the search process. After the DLR system detects the object, range-gate filtering can be performed accurately to significantly improve detection and discovery probabilities of the object. To further demonstrate a detailed feature about searching the object by the DLR system, a PD curve and its corresponding RD curve in a shaded block of FIG. 4A are locally enlarged. As shown in FIG. 4B, it can be seen that it cumulatively takes about 32 milliseconds for the object to enter the dwelling beam of the DLR system. During the 32 milliseconds, a minimum PD of the search of the object can reach 3.284 arc seconds, and an RD of the search at the same moment is 74.041 meters.

8) After the DLR system detects the object in the above constant elevation search process, it is assumed that one set of detection data {circumflex over (t)} and {circumflex over (ρ)} is obtained, where {circumflex over (t)} represents the detection moment, corresponding to 10:12:10.832 on Mar. 16, 2021 on Beijing time, and {circumflex over (ρ)} represents the detection range, which is equal to {circumflex over (ρ)}=1910478.253 meters and obtained based on the standard point of the object at the detection moment. The one set of detection data is detected by an 11^th chronological dwelling beam of the DLR system. An along-track error is determined based on the one set of detection data (a ranging delay involved in the calculation is given by the Hopfield model based on the ratio of 90%), to obtain an along-track error, namely, ϵ=−323.742 milliseconds. The along-track error is corrected for the initial quasi-mean elements of the object only based on the error ϵ generated by using the above one set of detection data. Tracking-specific guidance data of the object after moment {circumflex over (t)} during the current transit is generated based on a set of corrected initial quasi-mean elements, which is detection azimuths, detection elevations, and detection ranges that correspond to a series of detection moments. The detection azimuth and the detection elevation are used to determine a tracking direction of the DLR system, and the detection range is used to perform the range-gate filtering on tracking echo signals of the DLR system. The atmospheric refraction effect contained in the above tracking-specific guidance data is given by the Hopfield model based on the ratio of 90%. For any detection moment in the above tracking process, a PD of the tracking of the object at the moment is defined as an included angle between a tracking direction of the DLR system at the moment and a direction of the standard point of the object at the same moment, and an RD of the tracking of the object at the moment is defined as an absolute value of a difference between a detection range at the moment and a slant distance of the standard point of the object at the same moment. A situation of detecting the object by the DLR system in the above tracking process can be determined based on PDs of the tracking of the object at various detection moments. FIG. 5 is a coordinate graph with two Y-axes, which shows change curves of the PD of the tracking of the object and the RD of the tracking of the object over time in the above tracking process. It can be seen from FIG. 5 that for most periods in the entire tracking process, the PD of the tracking of the object is less than 10 arc seconds. The PD of the tracking slightly exceeds 10 arc seconds only in a short period at the beginning of the tracking process, and the RD of the tracking does not exceed 30 meters throughout the entire tracking process. This indicates that the DLR system can basically detect the object throughout the entire tracking process, and can perform the range-gate filtering accurately to greatly improve the detection and discovery probabilities of the object.

Based on the orbital parameters t_q, {right arrow over (r)}_q, {right arrow over ({dot over (r)})}_q, and ϵ_a and a high-precision perturbation calculation model, and taking into account the impact of the light travel time and the atmospheric refraction effect (given by the Hopfield model based on the ratio of 90%), the experiment generates orbital prediction data during the current transit, which are the detection azimuths, elevations, and ranges that correspond to a series of detection moments. For any of the above detection moments, a PD of the prediction of the object at the moment is defined as an included angle between a prediction direction of the object at the moment and a direction of the standard point of the object at the same moment, and an RD of the prediction of the object at the moment is defined as an absolute value of a difference between a detection range at the moment and a slant distance of the standard point of the object at the same moment. A situation of detecting the object by the DLR system under guidance of the above orbital predictions can be determined based on PDs of the prediction of the object at various detection moments. FIG. 6 is a coordinate graph with two Y-axes, which shows change curves of the PD of the prediction of the object and the RD of the prediction of the object over time during the current transit. From FIG. 6, it can be seen that the PD of the prediction of the object is greater than 10 arc seconds throughout the entire transit period and reaches tens of arc seconds in most of the transit period. The RD of the prediction of the object exceeds 100 meters in most of the transit period, and the PD of the prediction of the object exceeds 40 arc seconds in a period in which the RD of the prediction is less than 100 meters. This indicates that the DLR system is difficult to discover the object directly based on the guidance of the orbital predictions during the current transit, thus fully reflecting progressiveness of the technology in the present disclosure.

Embodiment 3

This embodiment selects SWARM B (NASA catalog number 39451) as an experimental object, which has an orbital height of 490 kilometers and is a laser calibration satellite. In the experiment, one set of orbital parameters of the object relative to the J2000 geocentric inertial system is generated based on detection data obtained by a precision tracking radar in China, namely t_q, {right arrow over (r)}_q, _q, and ϵ_a, as shown in Table 3:

TABLE 3
One set of orbital parameters of the SWARM B relative
to the J2000 geocentric inertial coordinate system
tq(Beijing time)	Feb. 28, 2021 10:35:25.752
{right arrow over (r)}q	xq	−176610.508	q	{dot over (x)}q	−6988.727
(m)	yq	223433.934	(m/s)	{dot over (y)}q	−2974.227
	zq	−6881097.319		żq	83.572
ϵa(m2/kg)	0.00027282

The experiment is based on the one set of orbital parameters and specified detection constraints of a DLR system, and a certain transit of the object is determined through high-precision perturbation calculation. The transit is a transit after 14 days starting from the epoch moment of the orbital parameters of the object. During the transit, search, and tracking experiments of the object are carried out.

1) A set of initial quasi-mean elements of the object at an intermediate moment during the current transit is obtained through perturbation propagation and conversion, as shown in Table 4:

TABLE 4
Initial quasi-mean elements of SWARM B at the
intermediate moment during the transit
T0 (Beijing time)	Mar. 14, 2021 14:56:08.900
ā0(m)	6875051.448	ι 0	87°.7608
Ω 0	198°.7996	ξ 0	−0.0000120676
η 0	−0.0012949970	λ 0	153°.2184

2) A moment of a theoretical closest approach point of the object during the current transit is calculated.

3) It is estimated that absolute value τ of a maximum along-track error of the initial quasi-mean elements is 30 seconds, and it is obtained that Δτ is 0.5 seconds. It is obtained through calculation that l=61. Then, apparent path cluster {Γ_k|61,−61} is generated.

4) After all undetectable apparent paths are removed, an apparent path cluster formed by the remaining apparent paths is still {Γ_k|61,−61}.

5) All apparent paths with unsatisfactory detectable arc lengths are removed, and an apparent path cluster formed by the remaining apparent paths is still {Γ_k|61,−61}.

6) Elevation scope [h_p, h_q] suitable for the DLR system to perform constant elevation search is determined, where

$h_p=5°·5001$ $h_q=81°·0686$

7) According to the processing methods in the steps 7, 8, and 9 in Embodiment 1, a total of three search elevations with non-intersecting effective search periods are obtained, and search-specific guidance data is generated for each search elevation.

FIG. 7 shows search elevations whose degrees are in an ascending order and their corresponding effective search periods in chronological order from left to right. As shown in FIG. 7, one of the search elevations is 28.0976 degrees, and a corresponding effective search period starts from 14:54:8.681 on Mar. 14, 2021 and ends at 14:55:8.378 on Mar. 14, 2021. A total of 307 sets of search-specific guidance data are generated for the elevation. The DLR system chronologically enables, based on each set of search-specific guidance data, the beam to dwell, to form a constant elevation search process.

FIG. 8A and FIG. 8B are coordinate graphs with two Y-axes, which respectively show changes to a PD of the search of the object and an RD of the search of the object over time in the above constant elevation search process. FIG. 8A shows change curves of the PD of the search and the RD of the search over time in the entire search process. From FIG. 8A, it can be seen that the object enters a dwelling beam of the DLR system during a very short period in the search process, and the RD of the search during this period is also very small, indicating that the DLR system can detect the object in the search process. After the DLR system discovers the object, range-gate filtering can be performed accurately to significantly improve detection and discovery probabilities of the object. To further demonstrate a detailed feature about detecting the object by the DLR system, a PD curve and its corresponding RD curve in a shaded block of FIG. 8A are locally enlarged. As shown in FIG. 8B, it can be seen that it cumulatively takes about 20 milliseconds for the object to enter the dwelling beam of the DLR system. During the 20 milliseconds, a minimum PD of the search of the object can reach 4.801 arc seconds, and an RD of the search at the same moment is 83.330 meters.

8) After the DLR system detects the object in the above constant elevation search process, it is assumed that one set of detection data {circumflex over (t)} and {circumflex over (ρ)} is obtained, where {circumflex over (t)} represents the detection moment, corresponding to 14:54:18.219 on Mar. 14, 2021 on Beijing time, and {circumflex over (ρ)} represents the detection range, which is equal to {circumflex over (ρ)}=955312.332 meters. The one set of detection data is detected by a 50^th chronological dwelling beam of the DLR system. Based on the one set of detection dat_q, an along-track error is determined to obtain an along-track error, namely Ε=20385.021 milliseconds. The along-track error is corrected for the initial quasi-mean elements of the object only based on the error ϵ generated by using the above one set of detection data. Tracking-specific guidance data of the object after moment {circumflex over (t)} during the current transit is generated based on a set of corrected initial quasi-mean elements.

FIG. 9 is a coordinate graph with two Y-axes, which shows change curves of the PD of the tracking of the object and the RD of the tracking of the object over time in the above tracking process. It can be seen from FIG. 9 that for most periods in the entire tracking process, the PD of the tracking of the object is less than 10 arc seconds. The PD of the tracking slightly exceeds 10 arc seconds only in a short period at the beginning of the tracking process, and the RD of the tracking of the object does not exceed 30 meters throughout the entire tracking process. This indicates that the DLR system can basically detect the object throughout the entire tracking process, and can perform the range-gate filtering accurately to greatly improve the detection and discovery probabilities of the object.

FIG. 10 is a coordinate graph with two Y-axes, which shows change curves of the PD of a prediction of the object and the RD of the prediction of the object over time during the current transit. From FIG. 10, it can be seen that the PD of the prediction of the object is greater than 1,000 arc seconds throughout the transit. The RD of the prediction of the object exceeds 5,000 meters in most periods, and the PD of the prediction of the object exceeds 12,000 arc seconds in a period in which the RD of the prediction is less than 5,000 meters. This indicates that it is difficult to discover the object directly based on guidance of the orbital predictions during DLR detection, thus fully reflecting progressiveness of the technology in the present disclosure.

Embodiment 4

This embodiment provides a computer-readable storage medium. The computer-readable storage medium stores a computer program, and the computer program enables a computer to perform the search and tracking method for full time-domain laser detection of space debris in Embodiment 1.

Embodiment 5

This embodiment provides an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor. The processor, when executing the computer program, implements the search and tracking method for full time-domain laser detection of space debris in Embodiment 1.

In the embodiments of the present disclosure, the computer storage medium may be a tangible medium that may contain or store a program used by or used in combination with an instruction execution system, apparatus or device. The computer storage medium may be, but is not limited to, an electronic, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus or device, or any combination thereof. More specific examples of the computer storage medium include: an electrical connection with one or more wires, a portable computer disk, a hard disk, a random access memory (RAM), a read-only memory (ROM), an erasable programmable ROM (an EPROM or a flash memory), an optical fiber, a portable compact disc read-only memory (CD-ROM), an optical storage device, a magnetic storage device, or any suitable combination thereof.

A person of ordinary skill in the art may be aware that, in combination with the examples described in the embodiments disclosed in the present disclosure, units and algorithm steps may be implemented by electronic hardware or a combination of computer software and electronic hardware. Whether the functions are performed by hardware or software depends on particular applications and design constraints of the technical solutions. A person skilled in the art may use different methods to implement the described functions for each particular application, but it should not be considered that the implementation goes beyond the scope of the present disclosure.

What is described above is merely the preferred implementations of the present disclosure, the scope of protection of the present disclosure is not limited to the above embodiments, and all technical solutions following the idea of the present disclosure fall within the scope of protection of the present disclosure. It should be noted that several modifications and adaptations made by those of ordinary skill in the art without departing from the principle of the present disclosure should fall within the scope of protection of the present disclosure.

Claims

1. A search and tracking method for full time-domain laser detection of space debris, implemented during each transit of a space debris object relative to a debris laser ranging (DLR) system, and comprising: based on a set of precision orbital parameters of an object that are determined based on data of precision radars and/or charge-coupled device (CCD) telescopes of a space monitoring network, and by taking an intermediate moment determined based on pre-obtained start and end moments of a current transit of the object as an initial time, generating a set of initial quasi-mean elements, and constructing an analytical perturbation model;discretizing an along-track error, obtaining, in combination with the analytical perturbation model, a series of new perturbation models for describing motion of all virtual objects, and forming an apparent path cluster by different apparent paths generated from the virtual objects during the current transit;determining a maximum detectable elevation hg and a minimum detectable elevation hp of each valid apparent path in the apparent path cluster, and then calculating search elevations of the DLR system;for each valid apparent path in the apparent path cluster, calculating a corresponding detection elevation based on a corresponding perturbation model and in consideration of an impact of the light travel time, and obtaining, by using a dichotomous approach based on a difference between the detection elevation and the search elevation, detection feature parameters of virtual objects corresponding to all apparent paths when the virtual objects rise to the search elevation;generating a series of sets of search-specific guidance data based on the detection feature parameters of the virtual objects, wherein each set of search-specific guidance data corresponds to one dwell of a laser beam emitted by the DLR system on the search elevation;performing, by the DLR system, multi-elevation search on the object based on the search-specific guidance data; andswitching to, by the DLR system, a tracking process after detecting the object, calculating an along-track error of the initial quasi-mean elements based on a few sets of detection data obtained before switching to the tracking process, to correct an along-track error of the initial quasi-mean elements, and generating tracking-specific guidance data based on a set of corrected initial quasi-mean elements to guide the DLR system to track the object.

2. The search and tracking method for full time-domain laser detection of space debris according to claim 1, further comprising: after constructing the analytical perturbation model, calculating, in combination with the constructed analytical perturbation model, a moment of a theoretical closest approach point of the object during the current transit; and for the apparent path cluster formed by the different apparent paths generated from the virtual objects during the current transit, determining visibility of each apparent path in the apparent path cluster based on the calculated moment of the theoretical closest approach point, removing all the invisible apparent paths and further removing all the apparent paths with an unsatisfactory detectable arc length, and taking the remaining apparent paths as valid apparent paths; andafter obtaining the detection feature parameters of the virtual objects corresponding to all the valid apparent paths when the virtual objects rise to the search elevation, performing continuation processing on a value of an azimuth in the detection feature parameters to make the value of the azimuth change continuously.

3. The search and tracking method for full time-domain laser detection of space debris according to claim 2, wherein the step of constructing the analytical perturbation model comprises: selecting the intermediate moment T0 of the current transit of the object as a reference moment, defining an orbital along-track error of the object during the current transit as an along-track error at the reference moment, and using the first-type nonsingular orbital elements as basic variables to construct an analytical perturbation model with the reference moment as the initial time:assuming that one set of known precision orbital parameters of the object are tq, {right arrow over (r)}q, {right arrow over ({dot over (r)})}q, and ϵa, where tq represents an epoch moment of the one set of orbital parameters, {right arrow over (r)}q and {right arrow over ({dot over (r)})}q respectively represent a position vector and a velocity vector of the object relative to an epoch geocentric inertial coordinate system, and ϵa represents an area-mass ratio of the object, performing perturbation propagation by using a numerical method and a precision dynamic model based on the known orbital parameters tq, {right arrow over (r)}q, {right arrow over ({dot over (r)})}q, and ϵa, performing propagation from the moment tq to the moment T0 to obtain a position vector {right arrow over (r)}0 and a velocity vector {right arrow over ({dot over (r)})}0 of the object relative to the epoch geocentric inertial coordinate system at the moment T0, performing conversion on {right arrow over (r)}0 and {right arrow over ({dot over (r)})}0 to obtain a set of initial quasi-mean elements of the object at the moment T0, and constructing, based on the initial quasi-mean elements, an analytical perturbation model with a following mathematical expression:assuming that α, i, Ω, ξ=e cos ω, η=−e sin ω, and λ=ω+M are first-type nonsingular osculating orbital elements of the object, where a represents the orbital semi-major axis, i represents the orbital inclination, Ω represents the orbital right ascension of ascending node, e represents the orbital eccentricity, ω represents an orbital argument of perigee, and M represents an orbital mean anomaly, obtaining following formulas:

4. The search and tracking method for full time-domain laser detection of space debris according to claim 3, wherein the step of forming the apparent path cluster by different apparent paths generated from the virtual objects during the current transit comprises: taking the along-track error of the initial quasi-mean elements as a random variable uniformly distributed within [−τ, τ], discretizing the random variable to convert a probability problem into a certainty problem, and defining a following formula:

5. The search and tracking method for full time-domain laser detection of space debris according to claim 4, wherein the step of obtaining the detection feature parameters of the virtual objects corresponding to all valid apparent paths when the virtual objects rise to the search elevation comprises: for one apparent path Γk in an apparent path cluster {Γk|m,n}, wherein {Γk|m,n} represents an apparent path cluster obtained by removing all the invisible apparent paths and further removing all the apparent paths with an unsatisfactory detectable arc length from {Γk|l,−l}, performing search calculation forwards from a moment {tilde over (t)}k of a closest approach point along a rising segment by using the dichotomous approach based on the model k, generating, in each step, a position vector of the object relative to a True Equator and Mean Equinox of Epoch (TEMEE) coordinate system at a corresponding moment, calculating a corresponding detection elevation based on the position vector and by considering the impact of the light travel time and an atmospheric refraction effect, determining a processing method for a next step based on a difference between the detection elevation and the search elevation, and finally obtaining detection feature parameters of the object when the object rises to the search elevation, comprising a detection moment, a detection range, and a detection azimuth; and repeating the above process for other apparent paths in the apparent path cluster to obtain detection feature parameters corresponding to all virtual objects.

6. The search and tracking method for full time-domain laser detection of space debris according to claim 5, wherein the step of generating the series of sets of search-specific guidance data comprises: for one apparent path Γk in the apparent path cluster {Γk|m,n}, wherein a detection moment, a detection range, and a detection azimuth of a corresponding virtual object when the virtual object rises to the search elevation are respectively tk, ρk, and A′k, where A′k represents an azimuth after the continuation processing, taking tk as a base point, and constructing, through cubic spline interpolation based on the known parameter ρk on the base point, an interpolation function in which a detection range ρ varies with a detection moment t; and taking A′k as a base point, and constructing, through the cubic spline interpolation based on the known parameter tk on the base point, an interpolation function in which a detection moment t varies with a detection azimuth A′, wherein the obtained interpolation functions are as follows:

7. The search and tracking method for full time-domain laser detection of space debris according to claim 1, wherein the step of performing, by the DLR system, multi-elevation search on the object based on the search-specific guidance data comprises: generating a series of gradually increasing search elevations according to a principle that search elevations are generated in an ascending order and effective search periods of the search elevations do not intersect with each other, generating search-specific guidance data for each search elevation, and starting, by the DLR system based on the search-specific guidance data for the corresponding search elevation, search from a start search elevation and gradually switching to a highest search elevation, whereinthe series of gradually increasing search elevations are generated according to following steps:1) selecting a start search elevation h* between a minimum searchable elevation hp and a maximum searchable elevation hq; presetting a constant χ, where χ>0 represents adjustment time required for the DLR system to switch from one search elevation to another search elevation; and setting k=1, hk=h*, and {tilde over (h)}=h*, wherein hk represents a to-be-generated search elevation;2) calculating an effective search period [{tilde over (t)}b, {tilde over (t)}e] of the DLR system on the elevation {tilde over (h)}, and setting hn={tilde over (h)} and hx=hq;3) setting h=(hn+hx)/2, calculating an effective search period [tb, te] of the DLR system on the elevation h, and performing calculation according to ζ=tb−{tilde over (t)}e−χ;4) if |ζ|<εt, setting k=k+1, hk=h, and {tilde over (h)}=h, and performing the step 2) to continue the calculation, wherein εt represents a given small value relative to χ; otherwise, performing next calculation;5) setting Δh=hq−h; and if ζ<0 and Δh<εh, ending the calculation process, wherein εh represents a given small value; otherwise, performing next calculation; and6) iif ζ<0, setting hn=h, and performing the step 3) to continue the calculation; otherwise, setting hx=h, and performing the step 3) to continue the calculation; andobtaining the series of gradually increasing search elevations hk according to the above setting method.

8. The search and tracking method for full time-domain laser detection of space debris according to claim 4, wherein the step of calculating the along-track error of the initial quasi-mean elements based on a few sets of detection data obtained before switching to the tracking process, to correct the along-track error of the initial quasi-mean elements, and generating the tracking-specific guidance data based on the corrected initial quasi-mean elements comprises: setting one set of detection data obtained by the DLR system before switching to the tracking process as {circumflex over (t)} and {circumflex over (ρ)}, which is obtained during a certain dwell of the laser beam on a certain search elevation, wherein {circumflex over (t)} represents a detection moment, {circumflex over (ρ)} represents a corresponding detection range;1) setting an iterative initial value ϵ=0 of the along-track error, and using an atmospheric correction model to calculate a ranging delay Δρ caused by tropospheric refraction;2) performing calculation according to λϵ=λ0+(n0+λ1)ϵ, replacing λ0 on the right side of the formula (6) with λϵ in the analytical perturbation model 0 to obtain a perturbation model ϵ, and obtaining a calculated value ρc of the detection range corresponding to the detection moment {circumflex over (t)} and a time change rate {dot over (ρ)}c of the calculated value based on the model ϵ and in consideration of an impact of the light travel time; and3) correcting and updating the along-track error according to ϵ=ϵ+({circumflex over (ρ)}−Δρ−ρc)/{dot over (ρ)}c, and then performing the step 2) to perform recalculation until along-track errors before and after the update meet certain convergence accuracy, to obtain one along-track error;for each of L sets of detection data obtained before the DLR system switches to the tracking process, wherein L≥1, generating one along-track error ϵ according to the above calculation process, wherein a finally determined along-track error is as follows:

9. A computer-readable storage medium, wherein the computer-readable storage medium stores a computer program, and the computer program enables a computer to perform the search and tracking method for full time-domain laser detection of space debris according to claim 1.

10. An electronic device, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor, when executing the computer program, implements the search and tracking method for full time-domain laser detection of space debris according to claim 1.

11. The computer-readable storage medium according to claim 9, wherein the search and tracking method for full time-domain laser detection of space debris further comprises: after constructing the analytical perturbation model, calculating, in combination with the constructed analytical perturbation model, a moment of a theoretical closest approach point of the object during the current transit; and for the apparent path cluster formed by the different apparent paths generated from the virtual objects during the current transit, determining visibility of each apparent path in the apparent path cluster based on the calculated moment of the theoretical closest approach point, removing all the invisible apparent paths and further removing all the apparent paths with an unsatisfactory detectable arc length, and taking the remaining apparent paths as valid apparent paths; andafter obtaining the detection feature parameters of the virtual objects corresponding to all the valid apparent paths when the virtual objects rise to the search elevation, performing continuation processing on a value of an azimuth in the detection feature parameters to make the value of the azimuth change continuously.

12. The computer-readable storage medium according to claim 11, wherein the step of constructing the analytical perturbation model comprises: selecting the intermediate moment T0 of the current transit of the object as a reference moment, defining an orbital along-track error of the object during the current transit as an along-track error at the reference moment, and using the first-type nonsingular orbital elements as basic variables to construct an analytical perturbation model with the reference moment as the initial time:assuming that one set of known precision orbital parameters of the object are tq, {right arrow over (r)}q, {right arrow over ({dot over (r)})}q, and ϵa, where tq represents an epoch moment of the one set of orbital parameters, {right arrow over (r)}q and {right arrow over ({dot over (r)})}q respectively represent a position vector and a velocity vector of the object relative to an epoch geocentric inertial coordinate system, and ϵa represents an area-mass ratio of the object, performing perturbation propagation by using a numerical method and a precision dynamic model based on the known orbital parameters tq, {right arrow over (r)}q, {right arrow over ({dot over (r)})}q, and ϵa, performing propagation from the moment tq to the moment T0 to obtain a position vector {right arrow over (r)}0 and a velocity vector {right arrow over ({dot over (r)})}0 of the object relative to the epoch geocentric inertial coordinate system at the moment T0, performing conversion on {right arrow over (r)}0 and {right arrow over ({dot over (r)})}0 to obtain a set of initial quasi-mean elements of the object at the moment T0, and constructing, based on the initial quasi-mean elements, an analytical perturbation model with a following mathematical expression:assuming that α, i, Ω, ξ=e cos ω, η=−e sin ω, and λ=ω+M are first-type nonsingular osculating orbital elements of the object, where a represents the orbital semi-major axis, i represents the orbital inclination, Ω represents the orbital right ascension of ascending node, e represents the orbital eccentricity, ω represents an orbital argument of perigee, and M represents an orbital mean anomaly, obtaining following formulas:

13. The computer-readable storage medium according to claim 12, wherein the step of forming the apparent path cluster by different apparent paths generated from the virtual objects during the current transit comprises: taking the along-track error of the initial quasi-mean elements as a random variable uniformly distributed within [−τ, τ], discretizing the random variable to convert a probability problem into a certainty problem, and defining a following formula:

14. The computer-readable storage medium according to claim 13, wherein the step of obtaining the detection feature parameters of the virtual objects corresponding to all valid apparent paths when the virtual objects rise to the search elevation comprises: for one apparent path Γk in an apparent path cluster {Γk|m,n}, wherein {Γk|m,n} represents an apparent path cluster obtained by removing the invisible apparent path and the apparent path with the unsatisfactory detectable arc length from {Γk|l,−l}, performing search calculation forwards from a moment {tilde over (t)}k of a closest approach point along a rising segment by using the dichotomous approach based on the model k, generating, in each step, a position vector of the object relative to a True Equator and Mean Equinox of Epoch (TEMEE) coordinate system at a corresponding moment, calculating a corresponding detection elevation based on the position vector and by considering the impact of the light travel time and an atmospheric refraction effect, determining a processing method for a next step based on a difference between the detection elevation and the search elevation, and finally obtaining detection feature parameters of the object when the object rises to the search elevation, comprising a detection moment, a detection range, and a detection azimuth; and repeating the above process for other apparent paths in the apparent path cluster to obtain detection feature parameters corresponding to all virtual objects.

15. The computer-readable storage medium according to claim 14, wherein the step of generating the series of sets of search-specific guidance data comprises: for one apparent path Γk in the apparent path cluster {Γk|m,n}, wherein a detection moment, a detection range, and a detection azimuth of a corresponding virtual object when the virtual object rises to the search elevation are respectively tk, ρk, and A′k, where A′k represents an azimuth after the continuation processing, taking tk as a base point, and constructing, through cubic spline interpolation based on the known parameter ρk on the base point, an interpolation function in which a detection range ρ varies with a detection moment t; and taking A′k as a base point, and constructing, through the cubic spline interpolation based on the known parameter tk on the base point, an interpolation function in which a detection moment t varies with a detection azimuth A′, wherein the obtained interpolation functions are as follows:

16. The computer-readable storage medium according to claim 9, wherein the step of performing, by the DLR system, multi-elevation search on the object based on the search-specific guidance data comprises: generating a series of gradually increasing search elevations according to a principle that search elevations are generated in an ascending order and effective search periods of the search elevations do not intersect with each other, generating search-specific guidance data for each search elevation, and starting, by the DLR system based on the search-specific guidance data for the corresponding search elevation, search from a start search elevation and gradually switching to a highest search elevation, whereinthe series of gradually increasing search elevations are generated according to following steps:1) selecting a start search elevation h* between a minimum searchable elevation hp and a maximum searchable elevation hq; presetting a constant χ, where χ>0 represents adjustment time required for the DLR system to switch from one search elevation to another search elevation; and setting k=1, hk=h*, and {tilde over (h)}=h*, wherein hk represents a to-be-generated search elevation;2) calculating an effective search period [{tilde over (t)}b, {tilde over (t)}e] of the DLR system on the elevation {tilde over (h)}, and setting hn={tilde over (h)} and hx=hq;3) setting h=(hn+hx)/2, calculating an effective search period [tb, te] of the DLR system on the elevation h, and performing calculation according to ζ=tb−{tilde over (t)}e−χ;4) if |ζ|<εt, setting k=k+1, hk=h, and {tilde over (h)}=h, and performing the step 2) to continue the calculation, wherein εt represents a given small value relative to χ; otherwise, performing next calculation;5) setting Δh=hq−h; and if ζ<0 and Δh<εh, ending the calculation process, wherein εh represents a given small value; otherwise, performing next calculation; and6) iif ζ<0, setting hn=h, and performing the step 3) to continue the calculation; otherwise, setting hx=h, and performing the step 3) to continue the calculation; andobtaining the series of gradually increasing search elevations hk according to the above setting method.obtaining the series of gradually increasing search elevations hk according to the above setting method.

17. The computer-readable storage medium according to claim 13, wherein the step of calculating the along-track error of the initial quasi-mean elements based on a few sets of detection data obtained before switching to the tracking process, to correct the along-track error of the initial quasi-mean elements, and generating the tracking-specific guidance data based on the corrected initial quasi-mean elements comprises: setting one set of detection data obtained by the DLR system before switching to the tracking process as {circumflex over (t)} and {circumflex over (ρ)}, which is obtained during a certain dwell of the laser beam on a certain search elevation, wherein {circumflex over (t)} represents a detection moment, {circumflex over (ρ)} represents a corresponding detection range;1) setting an iterative initial value ϵ=0 of the along-track error, and using an atmospheric correction model to calculate a ranging delay Δρ caused by tropospheric refraction;2) performing calculation according to λϵ=λ0+(n0+λ1)ϵ, replacing λ0 on the right side of the formula (6) with λϵ in the analytical perturbation model 0 to obtain a perturbation model ϵ, and obtaining a calculated value ρc of the detection range corresponding to the detection moment {circumflex over (t)} and a time change rate {dot over (ρ)}c of the calculated value based on the model ϵ and in consideration of an impact of the light travel time; and3) correcting and updating the along-track error according to ϵ=ϵ+({circumflex over (ρ)}−Δρ−ρc)/{dot over (ρ)}c, and then performing the step 2) to perform recalculation until along-track errors before and after the update meet certain convergence accuracy, to obtain one along-track error;for each of L sets of detection data obtained before the DLR system switches to the tracking process, wherein L≥1, generating one along-track error e according to the above calculation process, wherein a finally determined along-track error is as follows:

18. The electronic device according to claim 10, wherein the search and tracking method for full time-domain laser detection of space debris further comprises: after constructing the analytical perturbation model, calculating, in combination with the constructed analytical perturbation model, a moment of a theoretical closest approach point of the object during the current transit; and for the apparent path cluster formed by the different apparent paths generated from the virtual objects during the current transit, determining visibility of each apparent path in the apparent path cluster based on the calculated moment of the theoretical closest approach point, removing all the invisible apparent paths and further removing all the apparent paths with an unsatisfactory detectable arc length, and taking the remaining apparent paths as valid apparent paths; andafter obtaining the detection feature parameters of the virtual objects corresponding to all the valid apparent paths when the virtual objects rise to the search elevation, performing continuation processing on a value of an azimuth in the detection feature parameters to make the value of the azimuth change continuously.

19. The electronic device according to claim 18, wherein the step of constructing the analytical perturbation model comprises: selecting the intermediate moment T0 of the current transit of the object as a reference moment, defining an orbital along-track error of the object during the current transit as an along-track error at the reference moment, and using the first-type nonsingular orbital elements as basic variables to construct an analytical perturbation model with the reference moment as the initial time:assuming that one set of known precision orbital parameters of the object are tq, {right arrow over (r)}q, {right arrow over ({dot over (r)})}q, and ϵa, where tq represents an epoch moment of the one set of orbital parameters, {right arrow over (r)}q and {right arrow over ({dot over (r)})}q respectively represent a position vector and a velocity vector of the object relative to an epoch geocentric inertial coordinate system, and ϵa represents an area-mass ratio of the object, performing perturbation propagation by using a numerical method and a precision dynamic model based on the known orbital parameters tq, {right arrow over (r)}q, {right arrow over ({dot over (r)})}q, and ϵa, performing propagation from the moment tq to the moment T0 to obtain a position vector {right arrow over (r)}0 and a velocity vector {right arrow over ({dot over (r)})}0 of the object relative to the epoch geocentric inertial coordinate system at the moment T0, performing conversion on {right arrow over (r)}0 and {right arrow over ({dot over (r)})}0 to obtain a set of initial quasi-mean elements of the object at the moment T0, and constructing, based on the initial quasi-mean elements, an analytical perturbation model with a following mathematical expression:assuming that α, i, Ω, ξ=e cos ω, η=−e sin ω, and λ=ω+M are first-type nonsingular osculating orbital elements of the object, where a represents the orbital semi-major axis, i represents the orbital inclination, Ω represents the orbital right ascension of ascending node, e represents the orbital eccentricity, ω represents an orbital argument of perigee, and M represents an orbital mean anomaly, obtaining following formulas:

20. The electronic device according to claim 19, wherein the step of forming the apparent path cluster by different apparent paths generated from the virtual objects during the current transit comprises: taking the along-track error of the initial quasi-mean elements as a random variable uniformly distributed within [−τ, τ], discretizing the random variable to convert a probability problem into a certainty problem, and defining a following formula: