1. The celestial body target rotation estimation and three-dimensional reconstruction method is characterized by comprising the following steps of:

step one, modeling the motion of a celestial object:

under the celestial body inertial coordinate system O-XYZ with the celestial body center as the origin, the celestial body center coordinate is [0,0 ]]^TLet the coordinate of the radar be [ x ]₀(t)，y₀(t)，z₀(t)]^T；

The rotating shaft of the celestial body under an inertial coordinate system is taken as a unit vectorIn the direction of (a):

wherein a, b and c are unit vectors respectivelyProjection on three coordinate axes of a celestial body inertia coordinate system; y isThe included angle between the direction of the Z axis and the positive direction of the Z axis;is the angle between the X axis of the inertial coordinate system and OM, OM is the vectorProjection on the XOY plane of the inertial coordinate system; the p-th powder on the celestial bodyShoot point Q_pThe initial rectangular coordinate in the inertial coordinate system with the celestial body as the center is x_p0，y_p0，z_p0]^TP is 1,2, and P represents the number of scattering points; converting the inertial rectangular coordinate system with the celestial body center as the origin into a polar coordinate system to obtain the p-th scattering point Q_pThe corresponding initial polar coordinates may be expressed as [ alpha ]_p，β_p，H_p]^TWherein α is_pFor passing through the Z axis and Q of the celestial inertial coordinate system_pThe angle formed by the semi-plane of the upper part and the XOZ plane of the celestial body inertial coordinate system; beta is a_pIs OQ_pThe included angle between the Z axis and the positive direction of the Z axis; h_pThe distance from the p-th scattering point to the origin is the physical quantity to be searched; then, the rectangular coordinate under the celestial body inertial coordinate system corresponding to the pth scattering point is:

p scattering point around axis of rotationThe coordinates after rotation at the rotation speed ω are:

wherein r is [ a, b, ω ═ a, b, ω]^TThe rotation parameters of the celestial body are obtained; t is the slow time;

the slope course of the echo of the p-th scattering point after translational compensation is as follows:

wherein R is_T(t) is a compensated translation; in the formula (5), the parameters to be searched are a celestial body rotation parameter r and a height value H_p；

Step two: setting a parameter search space:

setting a search space of a celestial body rotation parameter r and a search space of scattering point heights;

step three: firstly, dividing a radar echo into a plurality of sub-apertures, obtaining a coherent accumulation result in each sub-aperture by utilizing GRFT, and then carrying out incoherent accumulation on results among the sub-apertures to obtain a final accumulation output value; wherein the accumulation process is represented as:

wherein s is_rRepresenting the echo after pulse compression; n is a radical of_cRepresenting the number of sub-apertures; t is_cRepresenting the time of each sub-aperture; λ is wavelength, c is speed of light, k is sub-aperture number;

then, in the three-dimensional parameter searching process, each group of rotation parameters r to be searched is fixed to be [ a, b, omega ]]^TFor each scattering point, the height value H_pP times of one-dimensional traversal search is carried out, a series of accumulated output values are obtained according to a formula (6) in each traversal search, and the maximum value of the P accumulated output values and the corresponding P height values are reserved; then, the P output maximum values are superposed to obtain an objective function s₁(r), expressed as:

finally, carrying out heuristic optimization on the rotation parameters in the parameter search space according to the formula (7), and obtaining the target function s of the formula (7)₁(r) the rotation parameter found when taking the maximum valueAnd height valueIs the true value of the target;

step four: using the obtained height values of P scattering pointsAnd reconstructing the three-dimensional shape of the celestial body.

2. The celestial object rotation estimation and three-dimensional reconstruction method of claim 1, wherein in step four, the three-dimensional reconstruction formula is:

wherein the content of the first and second substances,is the estimated value of the three-dimensional coordinates of the p-th scattering point.

3. The celestial object rotation estimation and three-dimensional reconstruction method of claim 1 or 2, wherein in step three, the non-coherent accumulation is performed by a Hybrid Generalized Radon Fourier Transform (HGRFT) based method.

4. The celestial object rotation estimation and three-dimensional reconstruction method of claim 1 or 2, wherein the heuristic optimization method is implemented by a particle swarm optimization algorithm.

5. The celestial object rotation estimation and three-dimensional reconstruction method of claim 1 or 2, wherein said rotation parameter search space is U ═ 1h,4h]×[-0.5，0.5]×[-0.5，0.5]Height search space of U_H＝[550m，800m]。

6. The celestial object rotation estimation and three-dimensional reconstruction method of claim 5, wherein the height search interval is 0.5 m.

Background

Inverse Synthetic Aperture Radar (ISAR) can image non-cooperative targets all-time and all-weather. The three-dimensional ISAR can reconstruct the three-dimensional structure of the target, obtain rich characteristic information of the target, and has a wide application prospect.

Three-dimensional ISAR imaging methods can be divided into two categories, namely an interferometric inverse synthetic aperture radar (InISAR) method and a single-antenna ISAR three-dimensional reconstruction method. Among other things, the InISAR method relies on a multi-base or multi-channel system that can form a synthetic aperture in the height dimension to achieve resolution in the height dimension. However, due to the influence of synchronization precision and coherence among different radars or channels, the method needs to invest higher hardware cost. The single-antenna ISAR three-dimensional reconstruction method obtains the synthetic aperture with the height dimension through the self movement of the target, the layout of an antenna array does not need to be considered, and the hardware complexity is low. Based on the ISAR image sequence of the single antenna, three-dimensional reconstruction can be realized by analyzing the projection equation of the target three-dimensional structure on the imaging plane, but the method is only limited to a cooperative target with known motion. In addition, based on the three-dimensional reconstruction theory of optical image processing, three-dimensional reconstruction can be realized according to the track matrix of the special display points in the ISAR image sequence, but the method needs to accurately calibrate the ISAR image in advance.

In recent years, deep space exploration has become a research hotspot. The practicability of the celestial body three-dimensional reconstruction method based on the InISAR is low due to the limitation of the number of ground-based radars. In addition, due to the long celestial body distance, the serious double-pass attenuation and the extremely low signal-to-noise ratio (SNR) of the echo, the ISAR image sequence is seriously interfered by noise, and the accurate calibration and three-dimensional reconstruction are difficult to realize. Therefore, the realization of the existing three-dimensional ISAR imaging method in celestial body observation has certain limitation.

Disclosure of Invention

In view of the above, the invention provides a method for estimating rotation and three-dimensional reconstruction of a celestial object based on radar long-time observation, which aims at the problem that three-dimensional ISAR imaging of the celestial object is difficult under low signal-to-noise ratio, and can realize joint estimation of a rotating shaft and a rotating speed of the celestial object under low signal-to-noise ratio and three-dimensional reconstruction of the celestial object.

1. The celestial body target rotation estimation and three-dimensional reconstruction method is characterized by comprising the following steps of:

step one, modeling the motion of a celestial object:

under the celestial body inertial coordinate system O-XYZ with the celestial body center as the origin, the celestial body center coordinate is [0,0 ]]^TLet the coordinate of the radar be [ x ]₀(t),y₀(t),z₀(t)]^T；

The rotating shaft of the celestial body under an inertial coordinate system is taken as a unit vectorIn the direction of (a):

wherein a, b and c are unit vectors respectivelyProjection on three coordinate axes of a celestial body inertia coordinate system; gamma isThe included angle between the direction of the Z axis and the positive direction of the Z axis;is the angle between the X axis of the inertial coordinate system and OM, OM is the vectorProjection on the XOY plane of the inertial coordinate system; p-th scattering point Q on celestial body_pThe initial rectangular coordinate in the inertial coordinate system with the celestial body as the center is x_p0,y_p0,z_p0]^TP is 1,2, …, and P represents the number of scattering points; converting the inertial rectangular coordinate system with the celestial body center as the origin into a polar coordinate system to obtain the p-th scattering point Q_pThe corresponding initial polar coordinates may be expressed as [ alpha ]_p,β_p,H_p]^TWherein α is_pFor passing through the Z axis and Q of the celestial inertial coordinate system_pThe angle formed by the semi-plane of the upper part and the XOZ plane of the celestial body inertial coordinate system; beta is a_pIs OQ_pThe included angle between the Z axis and the positive direction of the Z axis; h_pThe distance from the p-th scattering point to the origin is the physical quantity to be searched; then, the rectangular coordinate under the celestial body inertial coordinate system corresponding to the pth scattering point is:

p scattering point around axis of rotationThe coordinates after rotation at the rotation speed ω are:

wherein r is [ a, b, ω ═ a, b, ω]^TThe rotation parameters of the celestial body are obtained; t is the slow time;

the slope course of the echo of the p-th scattering point after translational compensation is as follows:

wherein R is_T(t) is a compensated translation; in the formula (5), the parameters to be searched are a celestial body rotation parameter r and a height value H_p；

Step two: setting a parameter search space:

setting a search space of a celestial body rotation parameter r and a search space of scattering point heights;

step three: firstly, dividing a radar echo into a plurality of sub-apertures, obtaining a coherent accumulation result in each sub-aperture by utilizing GRFT, and then carrying out incoherent accumulation on results among the sub-apertures to obtain a final accumulation output value; wherein the accumulation process is represented as:

wherein s is_rRepresenting the echo after pulse compression; n is a radical of_cRepresenting the number of sub-apertures; t is_cRepresenting the time of each sub-aperture; λ is wavelength, c is speed of light, k is sub-aperture number;

then, in the three-dimensional parameter searching process, each group of rotation parameters r to be searched is fixed to be [ a, b, omega ]]^TFor each scattering point, the height value H_pP times of one-dimensional traversal search is carried out, a series of accumulated output values are obtained according to a formula (6) in each traversal search, and the maximum value of the P accumulated output values and the corresponding P height values are reserved; then, the P output maximum values are superposed to obtain an objective function s₁(r), expressed as:

finally, carrying out heuristic optimization on the rotation parameters in the parameter search space according to the formula (7), and obtaining the target function s of the formula (7)₁(r) the rotation parameter found when taking the maximum valueAnd height valueIs the true value of the target;

step four: using the obtained height values of P scattering pointsAnd reconstructing the three-dimensional shape of the celestial body.

Preferably, in the fourth step, the formula of the three-dimensional reconstruction is as follows:

wherein the content of the first and second substances,is the estimated value of the three-dimensional coordinates of the p-th scattering point.

Preferably, in the third step, the non-coherent accumulation is implemented by using a hybrid generalized reed-solomon fourier transform (HGRFT) -based method.

Preferably, the heuristic optimization method is implemented by adopting a particle swarm optimization algorithm.

Preferably, the rotation parameter search space is U ═ 1h,4h]×[-0.5,0.5]×[-0.5,0.5]Height search space of U_H＝[550m,800m]。

Preferably, the height search interval is 0.5 m.

The invention has the following beneficial effects:

(1) the invention provides a target rotation parameter estimation and three-dimensional reconstruction method based on hybrid generalized Reed-Solomon Fourier transform (HGRFT). based on long-time accumulation, firstly, the spinning motion of a celestial body around a rotating shaft is modeled, and then coherent and non-coherent combined accumulation is realized by utilizing the HGRFT to the energy of an echo;

(2) when the motion model is matched with the real motion of the target, the accumulated peak value can be obtained, and the rotating shaft, the rotating speed and the height of a scattering point of the target can be quickly and accurately estimated according to the position of the peak value in the parameter searching space;

(3) the method can be combined with a heuristic optimization method to further improve the parameter searching efficiency.

Drawings

FIG. 1 shows a three-dimensional ISAR geometric model of a celestial object constructed in accordance with the present invention;

FIG. 2 shows a 1999KW4 star three-dimensional shape;

FIG. 3 shows the asteroid shape after meshing;

FIG. 4 shows the pulse compression results at a signal-to-noise ratio of-10 dB;

FIG. 5 shows the reconstructed asteroid shape;

fig. 6 shows the height value estimation deviation.

Detailed Description

The following describes in detail embodiments of the method of the present invention with reference to the accompanying drawings and examples.

If the radar transmits a chirp signal, the received echo can be expressed as:

wherein c is the speed of light; f. of_cIs the carrier frequency; λ c/f_cIs the wavelength; τ is the fast time; t is the slow time; t is_pIs the pulse width; k is a radical of_rIs the frequency modulation rate of the chirp signal; sigma_pThe amplitude of the p (p ═ 1,2, …) th scattering point echo; r_p(t) is the slope history of the p-th scattering point.

The echo after pulse compression can be represented as

Wherein, B is the signal bandwidth; a. the_p＝σ_pThe amplitude after pulse compression.

The implementation steps of the invention are as follows:

step one, modeling of target motion:

the distance between the target and the radar can be obtained according to the time delay of the received echo, and the azimuth angle and the pitch angle of the target relative to the radar can be determined according to the beam direction of the antenna. Then the celestial body center coordinate is [0,0 ] in the inertial coordinate system with the celestial body center as the origin]^TThe coordinate of the radar is [ x ]₀(t),₀(t),₀(t)]^TAs shown in fig. 1.

The rotating shaft of the celestial body under an inertial coordinate system is taken as a unit vectorIn the direction of (a) of (b), wherein a, b and c are unit vectors respectivelyProjection on three coordinate axes; gamma isThe included angle between the direction of the Z axis and the positive direction of the Z axis;is the angle between the X axis and OM, OM is the vectorProjection in the XOY plane. The initial rectangular coordinate of the P (P is 1,2, …, P) th scattering point on the celestial body under the inertial coordinate system with the celestial body as the center is [ x_p0,y_p0,z_p0]^TThen its corresponding initial polar coordinate can be expressed as [ alpha ]_p,β_p,H_p]^TWherein α is_pAnd beta_pThe coordinates of the p-th scattering point under the polar coordinate system need to be preset; h_pIs the height to be searched. Then the rectangular coordinate under the celestial body inertial coordinate system corresponding to the p-th scattering point is [ x ]_p0,y_p0,z_p0]^T

P scattering point around axis of rotationCoordinates after rotation at a rotation speed omega

Wherein r is [ a, b, ω ═ a, b, ω]^TIs the rotation parameter of the celestial body,

the slope course of the echo of the p-th scattering point after translational compensation is

Wherein R is_TAnd (t) is compensated translation, and can be obtained by the existing parameterized translation compensation method. In the formula (6), the parameters to be searched are a celestial body rotation parameter r and a three-dimensional height H_p(1,2,…,P)。

Step two: setting a parameter search space

Setting three-dimensional parameter search space U ═ omega_min,ω_max]×[a_min,a_max]×[b_min,b_max]For searching for a rotation parameter r, wherein (a)_min,b_min) And (a)_max,b_max) For search boundaries of celestial bodies' axes, omega_minAnd ω_maxA search boundary for the celestial body rotation speed; in addition, P identical one-dimensional parameter search spaces U are set_H＝[H_min,H_max]For searching the height of P scattering points, where H_minAnd H_maxIs a high search boundary.

Step three: based on hybrid generalized Reed-Solomon Fourier transform (HGRFT), searching parameter estimation values by combining traversal and heuristic optimization method, specifically:

the HGRFT method first divides the echo into N_cSub-apertures such that there is N within each sub-aperture_dA pulse, time T of each sub-aperture_c＝_dT_r. Wherein, mixing generalized raylsThe Fourier transform HGRFT method adopts the method of the patent with the application number of '202011196999.3', namely 'radar target detection method based on hybrid generalized Reed Fourier transform'. And obtaining coherent accumulation results by utilizing GRFT in each sub-aperture, and then carrying out incoherent accumulation on the results among the sub-apertures so as to obtain a final accumulation output value. In the present invention, the HGRFT-based accumulation process can be expressed as:

wherein k is the sub-aperture number.

In the invention, the parameter search is a process of three-dimensional search nested one-dimensional search. The heuristic optimization method can accelerate the multi-dimensional parameter searching process, so the heuristic optimization method is adopted in three-dimensional searching, and the traversal searching method is adopted in one-dimensional searching. The objective function of the search process is the accumulated output value of the HGRFT, and the specific steps are as follows.

Firstly, in the three-dimensional parameter searching process, each set of rotation parameters r to be searched is fixed to be [ a, b, omega ]]^TFor each scattering point, the height value H_pP times of one-dimensional traversal search is carried out, and P maximum output values of the HGRFTs and P corresponding height values are reserved. Next, the P maximum output values are superimposed and expressed as

And carrying out heuristic optimization on the rotation parameters in the three-dimensional parameter search space U according to the formula (8). When the searched rotation parameter and height value are searched to the real value of the target, the output of the equation (8) reaches the maximum value s₁(r) of (A). Final estimated values of the rotational parametersHeight value H from P scattering points_pCan be expressed as

Step four: reconstructing a three-dimensional shape of an object

After obtaining the height estimation values of P scattering pointsThe three-dimensional shape of the celestial body can be reconstructed into

Wherein the content of the first and second substances,is the estimated value of the three-dimensional coordinates of the p-th scattering point.

Example 1

In order to further verify the feasibility and the effectiveness of the invention, the rotation parameter estimation and three-dimensional reconstruction simulation of the proximal asteroid 1999KW4 primary star are carried out.

1999KW4 the main star has a rotation period of 2.8h and a diameter of 1.5km, and its three-dimensional shape is shown in FIG. 2. The surface of the small planet is subdivided at 200 m intervals, and the resulting three-dimensional shape consisting of 390 scattering points is shown in fig. 3. The simulation parameters are shown in Table 1

Table 1: simulation parameter table

Parameter(s)	Parameter value
		Carrier frequency	2380MHz
Bandwidth of	30MHz
		Sampling rate	30MHz
Pulse width	2us
		Pulse repetition frequency	1Hz
Number of subapertures	100
		Accumulation time	2000s
Radar coordinates	[1000m,1000m,500m]
		Rotating shaft	(0,0,1)

In the parameter searching process, a Particle Swarm Optimization (PSO) algorithm is adopted as a heuristic optimization method, wherein the number of PSO particles is set to be 100, the maximum iteration number is set to be 500, and the iteration threshold is set to be 10. The three-dimensional parameter search space is U ═ omega_min,ω_max]×[a_min,a_max]×[b_min,v_max]＝[1h,4h]×[-0.5,0.5]×[-0.5,0.5]Height search space of U_H＝[550m,800m]The height search interval is 0.5 m.

The noise was added to the echo after pulse compression to give a signal to noise ratio of-10 dB, and the result after pulse compression is shown in fig. 4, where it can be seen that the echo was completely drowned in noise. The method provided by the invention is adopted to search the rotation parameters, and the obtained estimated values are shown in the table 2. The reconstructed three-dimensional shape of the asteroid is shown in fig. 5, and the estimated deviations of the height values of 390 scattering points are shown in fig. 6, wherein 316 scattering points with estimated deviations smaller than 5 meters account for 81% of the total scattering points, 341 scattering points with estimated deviations smaller than 10 meters account for 87% of the total scattering points.

Table 2: estimation result of rotation parameter

	a	b	ω
				True value	0	0	2.8
Estimated value	0.0083315	0.005391	2.7952
				Deviation value	0.0083315	0.005391	0.0048

The invention can accurately estimate the rotation parameters and the three-dimensional shape of the celestial body, and the validity and the practicability of the invention are verified through simulation.