Patent ID: 12195209

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 inFIG.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 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 obtained detection data, 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 tq, {right arrow over (r)}q, {right arrow over ({dot over (r)})}q, and ϵa, where tqrepresents an epoch moment of the orbital parameters, {right arrow over (r)}qand {right arrow over ({dot over (r)})}qrespectively represent a position vector and a velocity vector of the object relative to an epoch geocentric inertial coordinate system, and ϵarepresents an area-mass ratio of the object (which is a ballistic coefficient related to atmospheric drag perturbation). In addition, Tband Teare 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:

T0=Tb+Te2

The first-type nonsingular orbital elements are used as basic variables, and an analytical perturbation model with T0as 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 tq, {right arrow over (r)}q, {right arrow over ({dot over (r)})}q, and ϵa. Propagation is performed from the moment tqto the moment T0to obtain position vector {right arrow over (r)}0and velocity vector {right arrow over ({dot over (r)})}0of the object relative to the epoch geocentric inertial coordinate system at the moment T0. Conversion is performed on {right arrow over (r)}0and {right arrow over ({dot over (r)})}0to obtain a set of initial quasi-mean elements of the object at the moment T0.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 T0, 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 a, 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:

a⁡(t)=a_0+as(1)(t)+as(2)(t)(1)i⁡(t)=ι_0+is(1)(t)+is(2)(t)(2)Ω⁡(t)=Ω_0+Ω1(t-T0)+Ωs(1)(t)+Ωs(2)(t)(3)ξ⁡(t)=ξ_0⁢cos[ω1(t-T0)]+η_0⁢sin[ω1(t-T0)]+ξs(1)(t)+ξs(2)(t)(4)η⁡(t)=η_0⁢cos[ω1(t-T0)]-ξ_0⁢sin[ω1(t-T0)]+ηs(1)(t)+ηs(2)(t)(5)λ⁡(t)=λ_0+(n_0+λ1)⁢(t-T0)+λs(1)(t)+λs(2)(t)(6)

In the above formulas, the left side represents osculating elements at moment t, ā0, ī0,Ω0,ξ0,η0, andλ0represent initial quasi-mean elements at the moment T0,n0=√{square root over (μ)}ā0−3/2represents an angular velocity of mean motion of the object, μ represents the geocentric gravitational constant, Ω1, ω1, and λ1represent coefficients of first-order long-term perturbation items of corresponding elements, as(1), is(1), Ωs(1), ξs(1), ηs(1), and λs(1)represent first-order short-period perturbation items of various elements, and as(2), is(2), Ωs(2), ξs(2), ηs(2), and λs(2)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*0of an approximate closest approach point is solved according to a following formula:

t0*=T0+λs-λ0λ.-u.s

In the above formula, λ0represents a value of λ(t) at the moment T0, which is the mean argument of latitude of the object at the moment T0and can be calculated according to the formula (6); {dot over (λ)}=n0+λ1represents a long-term change rate of λ(t); usand λsrespectively 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 T0; and {dot over (u)}srepresents a time change rate of us.

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

sin⁡(u~0-u0*)=(ξ~0⁢sin⁢u~0+η~0⁢cos⁢u~0)⁢sin2⁢θ(1-η~0⁢sin⁢u~0+ξ~0⁢cos⁢u~0)⁢cos⁢θ0(r/R-cos⁢θ)

In the above equation, {tilde over (ξ)}0and {tilde over (η)}0respectively represent values of ξ(t) and η(t) at the moment {tilde over (t)}0, 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)}0; θ0represents a latitude of the observation station with reference to orbit plane at the moment t*0; r represents a geocentric distance of the object at the moment {tilde over (t)}0; and R represents a geocentric distance of the observation station. Finally, the mean argument of latitude {tilde over (λ)}0of the object at the moment {tilde over (t)}0is calculated based on the elliptical motion relationship and ũ0, and the moment of the accurate closest approach point is calculated according to a following formula:

t~0=t0*+λ~0-λ0*λ.

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

Based on accuracy of the known orbital parameters tq, {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 T0is 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:

λ_k=λ_0+(n_0+λ1)⁢k⁢Δτ⁢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:

{Δτ=ετ≥2⁢εΔτ=τ/2τ<2⁢ε

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

l=int⁡(τ/Δτ)+1

In the constructed analytical perturbation model,λ0on 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 modelk. Each modelkdetermines one virtual object, and a virtual object corresponding to0is 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 Γkis uniquely determined by the modelk.

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 Γ0is 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 Γkand Γk+1are considered, a moment of an accurate closest approach point of Γkis 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, modelk+1is 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 Γ1, Γ2, and . . . , Γl.3) Two adjacent apparent paths Γkand Γk−1are considered, the moment of the accurate closest approach point of Γkis 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, modelk−1is 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 Γ−1, Γ−2, and . . . , Γ−l.

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

For one apparent path Γkin the apparent path cluster {Γk|l,−l}, moment {tilde over (t)}kof a closest approach point of the one apparent path is substituted into the modelkto 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 Γkin the apparent path cluster {Γk|i,j}, based on the modelkand 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)}kof 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 Γkin 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)}kof the closest approach point at an appropriate initial step by using a dichotomous approach based on the modelk. 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 hq. A maximum value is taken from the minimum detectable elevations of all the apparent paths, and set to hp. Then an elevation scope suitable for the DLR system to perform constant elevation search is determined based on closed interval [hp, hq], where hq≥hp. An available search elevation (ASE) of the DLR system can be obtained according to a following formula:

h=hp+β⁡(hq-hp)

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)}bof the target point relative to the Earth-fixed coordinate system and coordinates of the ground station. {right arrow over (r)}bcan be obtained through following coordinate conversion:

r→b=G⁡(t)⁢r→(7)

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)}bcan 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:

t=t′-ρ/c(8)

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)}bare interdependent, which are not known and can only be iteratively processed. During the iteration, initial condition t=t′ is set, {right arrow over (r)}bis 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 Γkin the apparent path cluster {Γk|m,n}, search calculation is performed forwards along the ascending segment from the moment {tilde over (t)}kof the closest approach point at an appropriate initial step by using the dichotomous approach based on the modelk. 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 Γkin the apparent path cluster {Γk|m,n}, an azimuth of a corresponding virtual object when the virtual object rises to the search elevation is Ak, where a value of Akis within interval [0,2π). The value of Akis redefined, and a defined value of Akis A′k.

k is considered to continuously change from m to n, and A′m=Amand ΔAk=Ak−Ak+1are set, where m−1≥k≥n. Subsequent A′kis calculated based on following three independent cases:1) If there is integer p that meets m−1≥p≥n, and ΔAp<−π, a following formula is obtained:

Ak′={Akk>pAk+2⁢πk≤p2) If there is integer p that meets m−1≥p≥n, and ΔAp>π, a following formula is obtained:

Ak′={Akk>pAk-2⁢πk≤p3) If there is no integer p meeting the above two cases, a following formula is obtained:

Ak′=Akm-1≥k≥n

A value of each A′kobtained above is no longer limited to the interval [0,2λ), but is consecutive while varying with k. A′kand Akcorrespond to a same physical orientation, and A′kcan be used to replace Akas the azimuth.

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

For one apparent path Γkin 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 tk, ρk, and A′k, where A′krepresents an azimuth after the continuation processing. k is considered to continuously change from m to n, and tkstrictly monotonically increases with k. Therefore, tkcan 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 ρkon the base point. A′kstrictly monotonically increases or decreases with k. Therefore, A′kcan 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 tkon the base point. The obtained interpolation functions are as follows:

ρ=ρ⁡(t)(9)t=t⁡(A′)(10)

A detection time scope when a real object rises to the search elevation is closed interval [tm, tn], 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:

Δψ=❘"\[LeftBracketingBar]"ψm-ψn❘"\[RightBracketingBar]"

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

Δ⁢A=❘"\[LeftBracketingBar]"Am′-An′❘"\[RightBracketingBar]"

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

d⁢ψ=cosh⁢dA∫ψmψnd⁢ψ=∫Am′An′cosh⁢dA

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

Δψ=cosh⁢Δ⁢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=int[Δψ(1-δ)⁢w]+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:

Δ⁢A*=Δ⁢A/N

A symbol factor is defined:

κ={1Am′<An′-1Am′>An′

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:

A~1*=Am′A~v+1*=A~v*+κ⁢Δ⁢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):

tv*=t⁡(A~v*)tv+1*=t⁡(A~v+1*)

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

A~v+1/2*=(A~v*+A~v+1*)/2

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

τv+1/2*=τ⁡(A~v+1/2*)

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

ρv+1/2*=ρ⁡(tv+1/2*)

The center azimuth Ã*v+1/2is 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:

Av+1/2*={A~v+1/2*-2⁢πA~v+1/2*≥2⁢πA~v+1/2*+2⁢πA~v+1/2*<0A~v+1/2*0≤A~v+1/2*<2⁢π

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/2is 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/2represents an azimuth of the direction of the beam, t*vand t*v+1respectively represent start and end dwell moments of the beam, and ρ*v+1/2represents 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*1=tm. 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=tn, and [tm, tn] 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 hpand maximum searchable elevation hqis 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, hk=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 hn={tilde over (h)} and hx=hqare set.3) It is set that h=(hn+hx)/2, effective search period [tb, te] of the DLR system on the elevation h is calculated through relevant calculation in the step 7, and calculation is performed according to ζ=tb−{tilde over (t)}e−χ.4) If |ζ|<εt, k=k+1, hk=h, and {tilde over (h)}=h are set, and the step 2) is performed to continue the calculation, where εtrepresents a given small value relative to χ; otherwise, next calculation is performed.5) It is set that Δh=hq−h. If ζ<0 and Δh<εh, the calculation process ends, where εhrepresents a given small value; otherwise, next calculation is performed.6) If ζ<0, hn=h is set, and the step 3) is performed to continue the calculation; otherwise, hx=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 hkare 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 data, 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λϵ=λ0+(n0+λ1)ϵ,λ0on the right side of the formula (6) is replaced withλϵand other formulas keep unchanged in the analytical perturbation model0to obtain perturbation modelϵ, and calculated value ρcof the detection range corresponding to the detection moment {circumflex over (t)} and time change rate {dot over (ρ)}cof 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 ϵ is generated according to the above calculation process. A finally determined along-track error is as follows:

ϵ_=∑ϵ/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λϵ=λ0+(n0+λ1)ϵ, replacingλ0on the right side of the formula (6) withλϵand keeping other formulas unchanged in the analytical perturbation model0to 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 tq, {right arrow over (r)}q,q, and ϵa, as shown in Table 1:

TABLE 1One set of orbital parameters of the AJISAI relativeto the J2000 geocentric inertial coordinate systemtq(Beijing time)Feb. 25, 2021 07:50:56.091{right arrow over (r)}qxq1059237.065q{dot over (x)}q−6942.578(m)yq4948397.066(m/s){dot over (y)}q1557.054zq6023858.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 T0of 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 tq, {right arrow over (r)}q, {right arrow over ({dot over (r)})}q, and ϵa. Propagation is performed from the moment tqto the moment T0to obtain position vector {right arrow over (r)}0and velocity vector {right arrow over ({dot over (r)})}0of the object relative to the J2000 geocentric inertial coordinate system at the moment T0. Further, {right arrow over (r)}0and {right arrow over ({dot over (r)})}0are converted into a set of initial quasi-mean elements of the object at the moment T0, as shown in Table 2.

TABLE 2Initial quasi-mean elements of AJISAI at the intermediatemoment during the transit 10:15:25.450T0(Beijing time)Mar. 16, 2021 10:15:25.450ā0(m)7866388.521ι050°.0079Ω0288°.1634ξ00.0010384276η0−0.0002674244λ032°.2032

Based on the initial quasi-mean elements, analytical perturbation model0describing 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 model0, two-step calculation is performed. In a first step, moment t*0of an approximate closest approach point is calculated based on T0. In a second step, moment {tilde over (t)}0of an accurate closest approach point is obtained through iterative calculation based on t*0. 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:

λ¯k=λ¯0+(n¯0+λ1)⁢k⁢Δ⁢τk=7,6,…,0,…,-6,-7

In the analytical perturbation model0,λkis used to replaceλ0in turn to obtain a series of perturbation modelsk. Each modelkdescribes 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 Γkin the apparent path cluster {Γk|7,−7}, based on the modelk, search calculation is carried out backwards along the descending segment from moment {tilde over (t)}kof 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 Γkin 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−3seconds. 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 hpis taken from the minimum detectable elevations of all the apparent paths, and minimum value hqis taken from the maximum detectable elevations of all the apparent paths, such that following values are obtained:

hp=5⁢°·5000hq=79⁢°·9082

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

7) Starting search elevation h*=15° is selected from [hp, hq]. 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 [{tilde over (t)}b, {tilde over (t)}e] for elevation {tilde over (h)} and effective search period [tb, te] 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−4degrees. 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.3shows search elevations whose degrees are in an ascending order and their corresponding effective search periods in chronological order from left to right. As shown inFIG.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 data, 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.4AandFIG.4Bare 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.4Ashows change curves of the PD of search and the RD of the search over time in the entire search process. FromFIG.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 ofFIG.4Aare locally enlarged. As shown inFIG.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 11thchronological 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.5is 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 fromFIG.5that 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 tq, {right arrow over (r)}q, {right arrow over ({dot over (r)})}q, and ϵaand 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.6is 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. FromFIG.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 tq, {right arrow over (r)}q,q, and ϵa, as shown in Table 3:

TABLE 3One set of orbital parameters of the SWARM B relativeto the J2000 geocentric inertial coordinate systemtq(Beijing time)Feb. 28, 2021 10:35:25.752{right arrow over (r)}qxq−176610.508q{dot over (x)}q−6988.727(m)yq223433.934(m/s){dot over (y)}q−2974.227zq−6881097.319żq83.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 4Initial quasi-mean elements of SWARM B at theintermediate moment during the transitT0(Beijing time)Mar. 14, 2021 14:56:08.900ā0(m)6875051.448ι087°.7608Ω0198°.7996ξ0−0.0000120676η0−0.0012949970λ0153°.21842) 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 [hp, hq] suitable for the DLR system to perform constant elevation search is determined, where

hp=5⁢°·5001hq=81⁢°·06867) 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.7shows search elevations whose degrees are in an ascending order and their corresponding effective search periods in chronological order from left to right. As shown inFIG.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.8AandFIG.8Bare 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.8Ashows change curves of the PD of the search and the RD of the search over time in the entire search process. FromFIG.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 ofFIG.8Aare locally enlarged. As shown inFIG.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 50thchronological dwelling beam of the DLR system. Based on the one set of detection data, 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.9is 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 fromFIG.9that 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.10is 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. FromFIG.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.