A method for analyzing the safety area of ​​a space uncontrolled tumbling target

Abstract

The invention discloses an analysis method for a safety area of ​​a spatial uncontrolled tumbling target, which comprises the following steps: step S1: measuring the angular velocity of the uncontrolled tumbling target and constructing an angular velocity coordinate system to obtain angular velocity point cloud data of the uncontrolled tumbling target in the angular velocity coordinate system; step S2: obtaining sparse nutation area point cloud data of the uncontrolled tumbling target; step S3: obtaining superimposed point cloud data, and performing homogenization processing on the superimposed point cloud data to obtain uniform nutation area point cloud data; step S4: performing mixed Gaussian modeling using the uniform nutation area point cloud data; step S5: obtaining a mixed Gaussian model of the uncontrolled tumbling target; and step S6: obtaining a safety area of ​​the uncontrolled tumbling target through the mixed Gaussian model. The invention takes into account the time-varying motion of the uncontrolled tumbling target, and describes the safety area of ​​the uncontrolled tumbling target through the mixed Gaussian model, thereby greatly improving the division accuracy of the safety area and facilitating the implementation of extremely close distance approach of the uncontrolled tumbling target.

Description

Technical Field

[0001] The invention belongs to the field of aerospace, and in particular relates to a method for analyzing a safe area of ​​a space uncontrolled tumbling target. Background Art

[0002] The rapid development of human space technology has spawned a large number of space missions. The number of space debris left in orbit is also increasing day by day, which seriously threatens the space environment and the safety of spacecraft in orbit. It is imperative to explore a method to systematically remove out-of-control targets in space orbit. Out-of-control space targets basically have no attitude adjustment ability and will roll under the influence of space perturbations. Therefore, they must be de-rotated before being captured. Existing contact or non-contact de-rotation methods all require maneuvering to a very close distance to the out-of-control space target to achieve a more ideal effect. In this process, it is crucial to avoid collision between the service satellite and the out-of-control space target.

[0003] In order to avoid collision with the uncontrolled tumbling target during the approach and derotation process, the safe area of ​​the uncontrolled tumbling target must be delineated. The commonly used methods for delineating the safe area are the ellipsoidal safe area and the safe corridor method. The ellipsoidal safe area is convenient for considering the influence of uncertainty and has a better effect in long-distance obstacle avoidance, but it cannot achieve the ideal obstacle avoidance effect in close-range maneuvers because it does not consider the complex shape of the uncontrolled tumbling target; the safe corridor method is suitable for docking tasks, but both the contact and non-contact derotation methods are inconvenient to implement in the conical safety corridor with the angular velocity of the uncontrolled tumbling target as the axis, and it is difficult to achieve the ideal derotation effect. Summary of the invention

[0004] The purpose of the present invention is to provide a method for analyzing the safety area of ​​a spatial uncontrolled rolling target, so as to solve the problem that when the uncontrolled rolling target approaches at a very close distance, the safety area does not take into account the time-varying motion of the uncontrolled rolling target and the division is not accurate enough.

[0005] The present invention adopts the following technical solution: a method for analyzing the safety area of ​​a space uncontrolled rolling target, which comprises the following steps:

[0006] Step S1: Measure the angular velocity of the uncontrolled rolling target and construct an angular velocity coordinate system. According to the original point cloud data of the uncontrolled rolling target in the body coordinate system and the measured angular velocity, a rotation matrix from the body coordinate system of the uncontrolled rolling target to the angular velocity coordinate system is established, and the angular velocity point cloud data of the uncontrolled rolling target in the angular velocity coordinate system is obtained.

[0007] Step S2: Rotate each point in the angular velocity point cloud data about the angular velocity axis by a random angle to obtain sparse nutation region point cloud data of the uncontrolled tumbling target.

[0008] Step S3: rotating the sparse nutation region point cloud data around the angular velocity axis and superimposing them to obtain superimposed point cloud data, and performing homogenization processing on the superimposed point cloud data to obtain uniform nutation region point cloud data.

[0009] Step S4: Use the uniform nutation region point cloud data to perform mixed Gaussian modeling, and use the log-likelihood function to evaluate the fitting effect.

[0010] Step S5: The lower bound of the log-likelihood function is continuously improved by the maximum expectation algorithm until the iteration termination condition is met, and the last iteration parameter value is taken as the model parameter to obtain the mixed Gaussian model of the uncontrolled rolling target.

[0011] Step S6: Obtain the safe area of ​​the uncontrolled rolling target through the mixed Gaussian model.

[0012] Furthermore, the deflection angle between the angular velocity point cloud data and the original point cloud data in step S1 is equal to the nutation angle of the uncontrolled tumbling target.

[0013] Furthermore, the probability density function of the mixed Gaussian model is:

[0014]

[0015] Where K is the number of selected Gaussian components, Θ = {Θ k =(π k ,μ k ,Σ k ), k = 1, 2, ... K} represents the parameter set of the mixed Gaussian model, π k is the weight of the kth Gaussian component, μ k is the mean of the kth Gaussian component, Σ k is the covariance of the kth Gaussian component, N(x|μ k ,Σ k ) indicates that the mean is μ k , the covariance is Σ k The probability density of the Gaussian distribution at coordinate point x.

[0016] Furthermore, the envelope of the safety area in step S6 is:

[0017] p(x|Θ)=p max

[0018] in,

[0019] Furthermore, the homogenization process in step S3 specifically consists of the following steps:

[0020] Step S301:

[0021] Sort the superimposed point cloud data from large to small according to the distance between each point and the origin.

[0022] Step S302:

[0023] Determine whether the distance between any two points and the origin satisfies R(j)-R(i)<c, where c is the expected point cloud spacing, R(j) represents the distance from the jth point to the origin, and R(i) represents the distance from the ith point to the origin.

[0024] Step S303:

[0025] When R(j)-R(i)≥c, the outer loop value is i+1, the inner loop value is j=i+1, and the process returns to step S302.

[0026] When R(j)-R(i)<c, determine whether the distance between any two points satisfies ||S(j)-S(i)||2<c, where S(i) represents the coordinates of the i-th point and S(j) represents the coordinates of the j-th point.

[0027] Step S304:

[0028] When ||S(j)-S(i)||2<c, delete point S(j) from the superimposed point cloud data, otherwise the inner loop value j+1 and return to step S302, and continue looping until the outer loop value i≥N, where N is the number of point clouds.

[0029] The beneficial effects of the present invention are:

[0030] 1. The present invention obtains the original point cloud data of the uncontrolled rolling target, and obtains the sparse nutation area point cloud data through processing, and then obtains the uniform nutation area point cloud data through homogenization processing, and then performs mixed Gaussian modeling on the uniform nutation area point cloud data, and determines the critical probability density value, then the equipotential surface of the obtained mixed Gaussian model can be used as the envelope surface of the safety area of ​​the uncontrolled rolling target;

[0031] 2. The present invention takes into account the time-varying motion of the uncontrolled rolling target and uses a mixed Gaussian model to characterize the safety area of ​​the uncontrolled rolling target, which greatly improves the accuracy of the division of the safety area and is conducive to the implementation of extremely close approach to the uncontrolled rolling target;

[0032] 3. The present invention takes the nutation area of ​​the uncontrolled tumbling target as the safety area, taking into account the time-varying motion of the uncontrolled tumbling target, and can ensure the safety of the de-rotation process when the nutation angle of the uncontrolled tumbling target changes slowly;

[0033] 4. The present invention uses the probability density equipotential surface of the mixed Gaussian model as the envelope surface of the safety area, and greatly improves the accuracy of the safety area division by utilizing the excellent fitting ability of the mixed Gaussian model, which can meet the requirements of tasks such as derotation of out-of-control tumbling targets in space at very close distances;

[0034] 5. The method proposed in the present invention has the advantages of the collision probability assessment method, and can compare the safety of each point according to the probability density value of the mixed Gaussian model, and is also convenient for defining the obstacle avoidance control force and planning the optimal trajectory;

[0035] 6. The homogenization algorithm of the present invention first classifies each point by its distance from the origin and then removes points that are too close, which reduces the amount of algorithm calculation and enables the nutation area to be updated in time when the nutation angle changes, thus meeting the real-time requirements of the task. BRIEF DESCRIPTION OF THE DRAWINGS

[0036] Figure 1 Schematic diagram of the original point cloud data of the present invention.

[0037] Figure 2 It is a schematic diagram of point cloud data of uniform nutation region of the present invention.

[0038] Figure 3 It is a schematic diagram of uniform nutation region point cloud data reconstructed by mixed Gaussian modeling of the present invention;

[0039] Figure 4 A schematic diagram of a rolling target safety area divided by the present invention;

[0040] Figure 5 Flow chart of the homogenization algorithm of the present invention. DETAILED DESCRIPTION

[0041] The present invention is described in detail below with reference to the accompanying drawings and specific embodiments.

[0042] The present invention discloses a method for analyzing a safe area of ​​a space uncontrolled rolling target, which comprises the following steps:

[0043] Step S1:

[0044] The angular velocity of the uncontrolled rolling target is measured and an angular velocity coordinate system is constructed. According to the original point cloud data of the uncontrolled rolling target in the body coordinate system and the measured angular velocity, a rotation matrix from the body coordinate system of the uncontrolled rolling target to the angular velocity coordinate system is established, and the angular velocity point cloud data of the uncontrolled rolling target in the angular velocity coordinate system is obtained.

[0045] Measuring the angular velocity ω of an uncontrolled rolling target in the Earth's equatorial inertial system I I =[ω x ,ω y ,ω z ] T , with angular velocity ω I The direction is the x-axis, that is, x H =ω I , at the vertical angular velocity ω ISelect any two orthogonal axes on the plane as the y and z axes, and let the vector corresponding to the selected y and z axes in the geocentric equatorial inertial system be y H =[y1,y2,y3] T 、z H =[z1,z2,z3] T , taking the target mass center as the origin O, construct the angular velocity coordinate system H = {O, x H ,y H ,z H}.

[0046] Before the implementation of the present invention, it is necessary to obtain the coordinate system of the uncontrolled rolling target body G = {O, x G ,y G ,z G}, such as Figure 1 shown.

[0047] According to the established angular velocity coordinate system, we first obtain the rotation matrix from the geocentric equatorial inertial system to the angular velocity coordinate system.

[0048]

[0049] Among them, [x I ,y I ,z I ]=E, E is the third-order unit matrix.

[0050] Then through the rotation matrix Get the vector of each coordinate axis of the target body coordinate system in the angular velocity coordinate system

[0051]

[0052]

[0053]

[0054] Then the rotation matrix from the target body coordinate system to the angular velocity coordinate system is for:

[0055]

[0056] Multiply each point in the original point cloud data P by the rotation matrix The main inertia axis of the original point cloud data deviates from its angular velocity axis, and the angular velocity point cloud data P' in the angular velocity coordinate system is obtained:

[0057]

[0058] Among them, P(i) is the coordinate vector of the i-th point in the original point cloud data P, and N0 is the number of points in the original point cloud data. According to the definition of nutation angle, it can be seen that the deflection angle between the original point cloud data and the angular velocity point cloud data in the angular velocity coordinate system is the nutation angle of the uncontrolled rolling target.

[0059] Step S2:

[0060] Each point in the angular velocity point cloud data is rotated around the angular velocity axis by a random angle to obtain sparse nutation region point cloud data of the uncontrolled tumbling target.

[0061] That is: rotate the angular velocity point cloud data in the angular velocity coordinate system around the x-axis of the angular velocity coordinate system by a random angle ψ i , and obtain the sparse nutation area point cloud data:

[0062] P”(i)=R(ψ i )P'(i)(i=1,2,…,N0)

[0063] in, P'(i) is the coordinate vector of the i-th point in the point cloud data P'.

[0064] Step S3:

[0065] The sparse nutation area point cloud data are rotated around the x-axis of the angular velocity coordinate system and superimposed to obtain superimposed point cloud data, and the superimposed point cloud data are homogenized to obtain uniform nutation area point cloud data, so that the point cloud density in the central area of ​​the uniform nutation area point cloud data is basically consistent with the point cloud density in the edge area.

[0066] The processed point cloud data P" is rotated several times around the x-axis of the angular velocity coordinate system and superimposed to obtain the superimposed point cloud data P*, and the superimposed point cloud data is homogenized to avoid excessive differences in the weights of each component when the mixed Gaussian model is used to model the data. Finally, the uniform nutation area point cloud data S* is obtained, as shown in Figure 2 shown.

[0067] The point cloud is homogenized, that is, points that are too close are deleted to ensure that the density of each part of the point cloud is basically the same. If the distance between each point is directly calculated to determine the points that are too close, it will generate a huge calculation load and cannot achieve a fast update of the nutation area. Therefore, the homogenization algorithm of the present invention classifies each point according to the distance from the origin, and only calculates the relative distance between points that are close to the origin. Figure 5 As shown, the main steps are:

[0068] (i) Sort the superimposed point cloud data P* from large to small according to the distance of each point from the origin to obtain the sorted point cloud data S, the distance R of each point from the origin, the number of point cloud data N, and the outer loop index value i is set to 1.

[0069] (ii) The inner loop index value j is set to j=i+1.

[0070] (III) Determine whether the condition R(j)-R(i)<c is satisfied, where c is the desired point cloud spacing, R(i) represents the distance from the i-th point to the origin, and R(j) represents the distance from the j-th point to the origin. When R(j)-R(i)≥c, the outer loop value is i+1, the inner loop value is j=i+1, and the process returns to step (III). If R(j)-R(i)<c is satisfied, determine whether the distance between the two points satisfies ||S(j)-S(i)||2<c, where S(i) represents the coordinates of the i-th point, and S(j) represents the coordinates of the j-th point. If ||S(j)-S(i)||2<c is satisfied, it is considered that the distance between the two points is too close, and the S(j) point is deleted from S, and the number of point cloud data N=N-1; if ||S(j)-S(i)||2<c is not satisfied, the S(j) point is retained. Determine whether the inner loop index value j satisfies j≥N. If so, the inner loop ends and the outer loop index value i=i+1. If not, the inner loop index value j=j+1 and this step is repeated. If R(j)-R(i)<c is not satisfied, the inner loop ends and the outer loop index value i=i+1.

[0071] (iv) Determine whether the outer loop index value i satisfies i≥N. If so, end the loop, homogenization is completed, and the remaining point cloud data is the uniform nutation area point cloud data S*; if not, return to step (ii) and continue the loop.

[0072] Step S4:

[0073] Mixed Gaussian modeling is performed using point cloud data from uniform nutation regions, and the log-likelihood equation is used to evaluate the fitting effect; mixed Gaussian modeling is a method of fitting any probability density function by the sum of a given number of Gaussian models. Due to its strong fitting ability and good function properties, it has been widely used in target modeling, uncertain propagation analysis and other fields, and has achieved certain results. The prior art models the target shape using a mixed Gaussian model, and uses the obtained mixed Gaussian model as a potential function to achieve obstacle avoidance. However, this method only uses the mixed Gaussian model as a means of evaluating the probability of collision, and cannot determine whether a collision occurs, so safety cannot be guaranteed. In addition, the spin of the target is not considered, and the requirements of the de-racemization task cannot be met.

[0074] Step S5:

[0075] The lower bound of the log-likelihood function is continuously improved through the maximum expectation algorithm until the iteration termination condition is met, and the last iteration parameter value is taken as the model parameter to obtain the mixed Gaussian model of the space out-of-control tumbling target.

[0076] Step S6:

[0077] The safe area of ​​the spatial uncontrolled tumbling target is obtained through the mixed Gaussian model, and the nutation area of ​​the uncontrolled tumbling target is taken as the safe area. The time-varying motion of the uncontrolled tumbling target is considered, and the safety of the de-rotation process can be guaranteed when the nutation angle of the uncontrolled tumbling target changes slowly.

[0078] Mixed Gaussian modeling is performed on the point cloud data S* of the uniform nutation region. According to experience, the number of Gaussian components is more appropriate to be 16-18. After determining the number of Gaussian components, the lower bound of the log-likelihood function of the model is continuously improved through the maximum expectation algorithm until the difference between the log-likelihood functions in two iterations is less than the termination error. Then the π obtained in the last iteration is selected. k ,μ k ,Σ k As the model parameter set Θ = {Θ k =(π k ,μ k ,Σ k ),k=1,2,…K}, the probability density function of the obtained mixed Gaussian model is:

[0079]

[0080] Where K is the number of selected Gaussian components, Θ = {Θ k =(π k ,μ k ,Σ k ), k = 1, 2, ... K} represents the parameter set of the mixed Gaussian model, π k is the weight of the kth Gaussian component, μ k is the mean of the kth Gaussian component, Σ k is the covariance of the kth Gaussian component, N(x|μ k ,Σ k ) indicates that the mean is μ k , the covariance is Σ k The probability density of the Gaussian distribution at coordinate point x.

[0081] The evaluation index is the logarithmic value of the joint distribution function of the uniform nutation region point cloud data S*, that is, the log-likelihood function:

[0082]

[0083] In the formula, x n represents the coordinates of the nth point in the point cloud, and L represents the likelihood value.

[0084] The larger the log-likelihood function value In(L), the better the fitting effect. Since there is still a summation term in the logarithmic term of the function, the maximum value cannot be obtained by differentiation, etc. Usually, according to Jensen inequality:

[0085] f[E(x)]≥E[f(x)]

[0086] This gives a lower bound on the log-likelihood function:

[0087]

[0088] In the formula, is the hidden parameter probability,

[0089] The lower bound Q(Θ) of the log-likelihood function is continuously improved by the maximum expectation algorithm. t ), until the iteration termination condition is met, the last iteration parameter value is taken as the model parameter, and the mixed Gaussian model of the point cloud data of the uniform nutation area is obtained. The nutation area reconstructed by the mixed Gaussian model is as follows Figure 3 As shown, one color represents the influence area of ​​a Gaussian component.

[0090] The main steps of the maximum expectation algorithm are:

[0091] (I) Iteration index value t = 1, input initial parameter Θ t .

[0092] (ii) According to parameter Θ t Calculate hidden parameter probabilities And give the lower bound Q(Θ) of the log-likelihood function t ).

[0093] (III) Find the solution that makes Q(Θ t ) takes the maximum value of the parameter value Θ t+1 .

[0094] (IV) Determine parameter value Θ t+1 Does it satisfy ||Θ t+1 -Θ t ||<α, where α is the iteration termination condition. If it is satisfied, the iteration is terminated and Θ is taken. t+1 As the parameter value of the mixed Gaussian model; if it is not satisfied, the index value t=t+1, and return to step (2) to continue iterating.

[0095] Determine the critical probability density value p max . Take the kth Gaussian component, whose mean is μ k , the Cholesky decomposition of the covariance matrix is ​​Σ k =GG T , any eigenvalue of G is δ i , the corresponding eigenvector is η i , take x0=μ i +3δ i (η i / ||η i|) as the critical probability density value p max , which can better envelop the point cloud data. The enveloping surface of the safe area is P(x|Θ)=p max , the safety constraint is P(x|Θ)<p max .

[0096] Critical probability density value p max If the safety zone is too small, it will be too large and the target of uncontrolled rolling cannot be approached at a very close distance; max If it is too large, a continuous envelope cannot be formed and collision avoidance cannot be achieved. Near the mean of any Gaussian component, the probability values ​​contributed by the other components are small, so it can be considered that:

[0097] p(μ k +σ|Θ)≈π k N(μ k +σ|μ k ,Σ k )

[0098] Where σ is a vector with a smaller bi-norm.

[0099] Covariance matrix Σ i is a symmetric positive definite matrix, and its Cholesky decomposition is Σ k =GG T , any eigenvalue of G is δ i , the corresponding eigenvector is η i , take x0=μ i +3δ i (η i / ||η i ||), then the probability density value of the Gaussian component at x0 is:

[0100]

[0101] According to the properties of the normal distribution, the confidence interval Ω of the normal distribution with confidence level γ satisfies:

[0102]

[0103] in:

[0104]

[0105] It represents the chi-square distribution value with n degrees of freedom and a significance level of 1-γ.

[0106] From the above formula, we can see that x0=μ i +3δ i (η i / ||η i||), the chi-square distribution value is 9. In this embodiment, the dimension of the point cloud data is 3, so the chi-square distribution degree of freedom n = 3. From the chi-square distribution table, it can be seen that the confidence γ is 0.9707 at this time, so x0 = μ i +3δ i (η i / ||η i ||) can ensure that more than 97% of the points in the point cloud corresponding to the selected Gaussian component are within the region Ψ, where:

[0107] Ψ={x:N(x|μ k ,Σ k )≥N(x0|μ k ,Σ k )}

[0108] The point cloud model has been homogenized, and the weights of the Gaussian components are approximately equal, so When P(x|Θ)≥p max With high confidence, it can be used as a safe area for uncontrolled tumbling targets in space, such as Figure 4 As shown, the safety constraint is P(x|Θ)<p max .

[0109] Measure and update the current angular velocity ω of the target. If in is the maximum permissible angular velocity value, return to step S1; if The de-rotation task is completed, the process ends, and the above steps are repeated to update the safe area until the de-rotation task is completed.

[0110] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims

1. A method for analyzing the safety area of ​​a spatial uncontrolled rolling target, characterized in that: It consists of the following steps: Step S1: Measure the angular velocity of the uncontrolled rolling target and construct an angular velocity coordinate system. According to the original point cloud data of the uncontrolled rolling target in the body coordinate system and the measured angular velocity, a rotation matrix from the body coordinate system of the uncontrolled rolling target to the angular velocity coordinate system is established, and the angular velocity point cloud data of the uncontrolled rolling target in the angular velocity coordinate system is obtained. Step S2: Rotate each point in the angular velocity point cloud data about the angular velocity axis by a random angle to obtain sparse nutation region point cloud data of the uncontrolled tumbling target. Step S3: rotating the sparse nutation region point cloud data around the angular velocity axis and superimposing them to obtain superimposed point cloud data, and performing homogenization processing on the superimposed point cloud data to obtain uniform nutation region point cloud data. Step S4: Use the uniform nutation region point cloud data to perform mixed Gaussian modeling, and use the log-likelihood function to evaluate the fitting effect. Step S5: The lower bound of the log-likelihood function is continuously improved by the maximum expectation algorithm until the iteration termination condition is met, and the last iteration parameter value is taken as the model parameter to obtain the mixed Gaussian model of the uncontrolled rolling target. Step S6: Obtain the safe area of ​​the uncontrolled rolling target through the mixed Gaussian model.

2. The method for analyzing the safety area of ​​a spatial uncontrolled rolling target according to claim 1 is characterized in that: The deflection angle between the angular velocity point cloud data and the original point cloud data in step S1 is equal to the nutation angle of the uncontrolled tumbling target.

3. The method for analyzing the safety area of ​​a spatial uncontrolled rolling target according to claim 2 is characterized in that: The probability density function of the mixed Gaussian model is: Where K is the number of selected Gaussian components, Θ = {Θ k =(π k ,μ k ,Σ k ), k=1,2,…K} represents the parameter set of the mixed Gaussian model, π k is the weight of the kth Gaussian component, μ k is the mean of the kth Gaussian component, Σ k is the covariance of the kth Gaussian component, N(x|μ k ,Σ k ) indicates that the mean is μ k , the covariance is Σ k The probability density of the Gaussian distribution at coordinate point x.

4. The method for analyzing the safety area of ​​a spatial uncontrolled rolling target according to claim 3 is characterized in that: The envelope of the safety area in step S6 is: p(x|Θ) p max in, 5. A method for analyzing the safety area of ​​a spatial uncontrolled rolling target according to any one of claims 1 to 4, characterized in that: The homogenization process in step S3 specifically consists of the following steps: Step S301: Sort the superimposed point cloud data from large to small according to the distance between each point and the origin. Step S302: Determine whether the distance between any two points and the origin satisfies R(j)-R(i)<c, where c is the expected point cloud spacing, R(j) represents the distance from the jth point to the origin, and R(i) represents the distance from the ith point to the origin. Step S303: When R(j)-R(i)≥c, the outer loop value is i+1, the inner loop value is j=i+1, and the process returns to step S302. When R(j)-R(i)<c, determine whether the distance between any two points satisfies ||S(j)-S(i)||2<c, where S(i) represents the coordinates of the i-th point and S(j) represents the coordinates of the j-th point. Step S304: When ||S(j)-S(i)||2<c, delete point S(j) from the superimposed point cloud data, otherwise the inner loop value j+1 and return to step S302, and continue looping until the outer loop value i≥N, where N is the number of point clouds.

Images