Radar attitude pointing method under large double-base angle system

Abstract

The invention discloses a radar attitude pointing method under a large double-base angle system, and the method comprises the steps: inputting a map image, and automatically recognizing a designated region according to a color threshold value; inputting six satellite orbit elements, UTC time and a target position; calculating the position of the satellite in a WGS-84 coordinate system; calculating a linear beam vector from the satellite to the target; determining three axes under a satellite body coordinate system; calculating an attitude conversion matrix according to a maximum energy principle, and calculating an attitude Euler angle from the attitude conversion matrix; and calculating the beam footprint overlapping ratio. According to the method, the vacancy of the attitude control method of the satellite-borne bistatic radar under the large bistatic angle system is filled, and the radar performance is improved; meanwhile, due to the fact that the observation angle difference is large, multi-angle target information can be conveniently obtained, the dimension and richness of data are enhanced, and higher-quality target analysis and evaluation are provided.

Description

TECHNICAL FIELD

[0001] The present application relates to a radar attitude pointing method, in particular to a radar attitude pointing method under a large bistatic angle mode. BACKGROUND

[0002] A spaceborne bistatic radar system is configured with two or more radar stations, which respectively transmit or receive radar waves from different angles to achieve synergistic effect.

[0003] The synchronization of spatial bistatic beams is a key problem. The purpose of spatial synchronization is to ensure that the beams of the active transmitting and passive receiving two radars, which are geometrically placed separately, are directed to the same target ground during radar operation, and have sufficient "footprint" overlap on the ground to ensure the quality of the receiving radar imaging. In the application of InSAR, Lou Liangsheng et al. proposed spatial synchronization methods of maximum coherence and maximum energy, as shown in Figure 1 The synchronization scheme of the maximum energy method can make the auxiliary star beam footprint coincide with the main star beam footprint by biasing the auxiliary star guiding law, without considering the beam pointing error, so that the auxiliary star beam center is directed to the center of the mapping band, and the auxiliary star received echo energy is maximized. The maximum coherence method does not require the main and auxiliary star beam centers to be directed to the same point on the ground at the same time, but requires the main and auxiliary stars to respectively irradiate the beam footprint in the same mapping band according to their respective signal and pitch two-dimensional guiding laws, and the synchronized beams do not completely overlap in the azimuth direction.

[0004] In the application of space-based bistatic radar observation of targets, compared with the distance between the two stars in InSAR, the interstellar distance is sometimes far, and there will be a large bistatic angle, but how to specify the attitude when observing in the large bistatic angle is an important problem. SUMMARY

[0005] In view of the large bistatic angle mode formed by the spaceborne bistatic radar during observation, the present application proposes a radar attitude pointing method under a large bistatic angle mode, which reduces the influence of single direction interference and improves the stability and reliability of the system in a complex electronic countermeasure environment. At the same time, using data of different angles, more comprehensive information about the target can be obtained, which helps to improve the accuracy of target recognition.

[0006] The purpose of the present application is achieved by the following technical solutions.

[0007] A radar attitude pointing method under a large bistatic angle mode, the steps comprising:

[0008] 1) Input the map image, and automatically identify the specified area according to the color threshold value;

[0009] 2) Input the satellite orbit six elements, UTC time and target position;

[0010] 3) Calculate the position of the satellite in WGS-84 coordinate system;

[0011] 4) Calculate the line-of-sight vector of the satellite to the target;

[0012] 5) Determine the three axes in the satellite body coordinate system;

[0013] 6) Calculate the attitude transformation matrix according to the principle of maximum energy, and calculate the attitude Euler angles from the attitude matrix;

[0014] 7) Calculate the beam footprint overlap rate.

[0015] Step 1) is specifically: 1-1) Read and convert the image: convert from RGB color space to HSV color space;

[0016] 1-2) Color threshold setting: set threshold in HSV color space to identify specific areas;

[0017] 1-3) Apply color filtering: create a mask based on the set threshold to isolate the river part;

[0018] 1-4) Extract target area coordinates: apply the mask to the original image to extract the coordinates of the river area.

[0019] Step 2) is specifically: 2-1) Use orbital six elements to describe the orbital coordinate system, orbital six elements are semi-major axis ( ), eccentricity ( ), inclination ( ), right ascension of the ascending node ( ), argument of periapsis ( ), true anomaly ( );

[0020] 2-2) When the satellite is in motion, the time t passed needs to be input to calculate the new true anomaly ( );

[0021] Step 2-2) is specifically: 2-2-1) Calculate the initial argument of periapsis :

[0022] Given the original true anomaly , calculate the argument of periapsis :

[0023] ;

[0024] 2-2-2) Calculate the mean anomaly :

[0025] ;

[0026] 2-2-3) Calculation of the orbital period :

[0027] Orbital period is the time it takes for the satellite to complete one revolution around the central body and is calculated by Kepler's third law:

[0028] ;

[0029] where is the standard gravitational parameter of the central body;

[0030] 2-2-4) Update the mean anomaly :

[0031] ;

[0032] As time increases, the mean anomaly will increase at a fixed rate, which is determined by the orbital period ;

[0033] 2-2-5) Solve Kepler's equation to obtain the eccentric anomaly using iterative approximation ;

[0034] ;

[0035] 2-2-6) Update the true anomaly ;

[0036] .

[0037] Step 3) Specifically: 3-1) Calculate the position of the satellite in the inertial coordinate system from the orbital elements , the calculation method is as follows:

[0038] Orbital plane determination:

[0039] ;

[0040] ;

[0041] Calculate the coordinates according to the rotation matrix:

[0042] ;

[0043] ;

[0044] ;

[0045] ;

[0046] 3-2) Conversion of inertial coordinate system to WGS-84 coordinate system:

[0047] According to the UTC time, the angle is calculated using the Greenwich Mean Sidereal Time, and the rotation matrix R is used to convert the WGS-84 coordinate system coordinates to the inertial coordinate system coordinates, and the calculation method is as follows:

[0048] The rotation angle is calculated according to the Greenwich Mean Sidereal Time calculation formula

[0049]

[0050] Where JD is the input UTC time corresponding to the Julian day, J2000 corresponds to the Julian day 2451545, and the last two terms are the precession correction terms;

[0051]

[0052]

[0053] Step 4) is specifically:

[0054] 4-1) Define beam direction: in the satellite body coordinate system, assume that the beam is emitted along the Z axis, that is, the beam direction in the body coordinate system is represented as

[0055] 4-2) Coordinate conversion: use the satellite attitude to convert the beam direction from the WGS-84 coordinate system to the orbit coordinate system VVLH;

[0056] The inverse matrix of the coordinate axis is:

[0057]

[0058] The transformation matrix around x, y, z is as follows:

[0059]

[0060]

[0061]

[0062]

[0063]

[0064] Step 5) is specifically:

[0065] 5-1) Determine the double-star axis pointing:

[0066] ​​​​​​​​​​​The radar beam center always points to the ground target particle, that is, the beam center vector of the bistatic radar in the ground-fixed system at any time , is expressed as

[0067] ;

[0068] ;

[0069] wherein, represents the position of satellite 1, represents the position of satellite 2, and the positions of the satellites in the ground-fixed system can be calculated through six orbital parameters and UTC time, represents the target position;

[0070] The beam center vector of the bistatic radar in the ground-fixed system at any time is converted to the orbital coordinate system and unitized to obtain , , , :

[0071] is the inverse matrix of the coordinate axis:

[0072] ;

[0073] The transformation matrix around x, y, and z is as follows:

[0074] ;

[0075] ;

[0076] , ;

[0077] 5-2) Determining the pointing direction of the double-satellite axis:

[0078] Due to the principle of maximum energy, the footprints of the double-satellite beams need to be coincident, and therefore it is set that the double-satellite axis points to the direction of the running speed of the middle position of the double satellite;

[0079] 5-3) Determining the pointing direction of the double-satellite axis:

[0080] After the axis and the axis are determined, the unit vector of the axis at time t, and can be obtained according to the right-hand rule.

[0081] Step 6) specifically involves calculating the Euler angles of the attitude generated by the transformed reference coordinate system coordinates and the original body coordinate system along the z-axis:

[0082] ,in, This is called the direction cosine matrix;

[0083] ;

[0084] For the direction cosine matrix ZYX Euler angles can be resolved as follows:

[0085] ;

[0086] ;

[0087] ;

[0088] Among them, here It is a matrix The line, number Column elements;

[0089] Step 7) specifically refers to:

[0090]

[0091] in, The overlapping area of ​​the two beam footprints. The area of ​​the first beam footprint, Let be the area of ​​the second beam footprint.

[0092] Compared with existing technologies, the advantages of this invention are: this invention fills the gap in attitude control methods for spaceborne bistatic radar under large bistatic angle conditions, improving radar performance; at the same time, due to the large difference in observation angles, it is easy to acquire multi-angle target information, enhance the dimensionality and richness of the data, and provide higher quality target analysis and evaluation. Attached Figure Description

[0093] Figure 1 This is a schematic diagram of the maximum energy method and the maximum coherence method.

[0094] Figure 2 This is the main flowchart of the present invention.

[0095] Figure 3 This is a schematic diagram of beam footprint overlap in this invention.

[0096] Figure 4 Set the track parameters table.

[0097] Figure 5The edge extraction result image.

[0098] Figure 6 The attitude result change curve diagram obtained by using matlab simulation.

[0099] Figure 7 The STK beam simulation result diagram. DETAILED DESCRIPTION

[0100] The application will be described in detail below with reference to the accompanying drawings and specific embodiments.

[0101] As Figure 2 shown, a radar attitude pointing method under a large dual-base angle system includes the following steps:

[0102] (1) Input the map image, and automatically identify the specified area according to the color threshold;

[0103] (1-1) Read and convert the image: convert from the RGB color space to the HSV color space; H (hue) in the HSV color space can more easily distinguish colors, while S (saturation) and V (brightness) help identify the intensity and lightness of the color.

[0104] (1-2) Color threshold setting: set the threshold in the HSV color space to identify a specific area. For example, a river may appear in a specific blue or dark color, which can be selected by adjusting the threshold.

[0105] (1-3) Apply color filtering: create a mask based on the set threshold to isolate the river part.

[0106] (1-4) Extract the target area coordinates: apply the mask to the original image to extract the coordinates of the river area.

[0107] (2) Input the orbital elements of two satellites, UTC time, and target position;

[0108] (2-1) The orbital elements are usually used to describe the orbital coordinate system, and the orbital elements are semi-major axis ( ), eccentricity ( ), inclination ( ), right ascension of the ascending node ( ), argument of periapsis ( ), and true anomaly ( ).

[0109] (2-2) When the satellite is in motion, the time t passed needs to be input to calculate the new true anomaly ( );

[0110] (2-2-1) Calculate the initial argument of periapsis :

[0111] Known true anomaly , calculate the eccentric anomaly :

[0112]

[0113] (2-2-2) Calculate the mean anomaly :

[0114]

[0115] (2-2-3) Calculate the orbital period :

[0116] Orbital period is the time required for a satellite to complete one revolution around the central celestial body, which can be calculated by Kepler's third law:

[0117]

[0118] where is the standard gravitational parameter of the central celestial body.

[0119] (2-2-4) Update the mean anomaly :

[0120]

[0121] This means that as time increases, the mean anomaly will increase at a fixed rate, and this rate is determined by the orbital period .

[0122] (2-2-5) Solve the Kepler equation to obtain the eccentric anomaly by using iterative approximation ;

[0123]

[0124] (2-2-6) Update the true anomaly ;

[0125]

[0126] (3) Calculate the position of the satellite in the WGS-84 coordinate system

[0127] (3-1) Calculate the position of the satellite in the inertial coordinate system from the orbital elements , the calculation method is as follows:

[0128] Orbital plane determination:

[0129]

[0130]

[0131] According to the rotation matrix, the coordinates are calculated:

[0132]

[0133]

[0134]

[0135]

[0136] (3-2) Conversion of inertial coordinate system to WGS-84 coordinate system

[0137] According to the UTC time, the angle is calculated using the Greenwich Mean Sidereal Time, and the coordinates of the WGS-84 coordinate system are converted to the coordinates of the WGS-84 coordinate system by the rotation matrix R, and the calculation method is as follows:

[0138] The rotation angle is calculated according to the Greenwich Mean Sidereal Time calculation formula

[0139]

[0140] where JD is the input UTC time corresponding to the Julian day, J2000 corresponds to the Julian day 2451545, and the last two terms are the precession correction terms.

[0141]

[0142]

[0143] (4) Calculate the satellite-to-target line-of-sight vector;

[0144] (4-1) Define the beam direction: in the satellite body coordinate system, assume that the beam is emitted along the Z axis, i.e. the beam direction in the body coordinate system can be represented as .

[0145] (4-2) Coordinate conversion: use the satellite attitude to convert the beam direction from the WGS-84 coordinate system to the orbit coordinate system VVLH.

[0146] The inverse matrix of the coordinate axes is:

[0147]

[0148] The transformation matrix around x, y, z is as follows:

[0149] ​​

[0150]

[0151]

[0152]

[0153]

[0154] (5) Determine the three-axis pointing of satellite body coordinate system

[0155] (5-1) Determine the two-axis pointing of satellite

[0156] The radar beam center always points to the ground target particle, that is, the beam center vector of the bistatic radar at any time in the earth-fixed system 、 can be expressed as

[0157]

[0158]

[0159] wherein, represents the position of satellite 1, represents the position of satellite 2, and the positions of satellites in the earth-fixed system can be calculated by the orbital six elements and UTC time, represents the target position.

[0160] Convert the beam center vector of the bistatic radar at any time in the earth-fixed system 、 to the orbital coordinate system and unitize to obtain 、 :

[0161] The inverse matrix of the coordinate axis is:

[0162]

[0163] The transformation matrix around x, y, z is as follows:

[0164]

[0165]

[0166] ,

[0167] (5-2) Determine the two-axis pointing of satellite

[0168] ​​Due to the principle of maximum energy, the dual-beam footprints need to be coincident, so the dual-beam The axes are all pointing to the direction of the velocity of the dual-beam.

[0169] (5-3) Determine the dual-beam The axes are pointing to

[0170] After determining the axes of the dual-beam and the axes of the dual-beam At time t, the unit vector of the axis of the dual-beam The unit vector of the axis of the dual-beam and the unit vector of the axis of the dual-beam .

[0171] (6) Attitude conversion matrix, calculate the Euler angle from the attitude matrix.

[0172] The Euler angle of the attitude generated by the converted reference coordinate system coordinates and the original body coordinate system z-axis direction

[0173] From the above, where, is called the direction cosine matrix.

[0174]

[0175] For the direction cosine matrix , the ZYX Euler angle can be analyzed in the following way:

[0176]

[0177]

[0178]

[0179] where, here is the element of the matrix in the i-th row and j-th column. (7) Calculate the beam footprint overlap rate;

[0180]

[0181] where, is the overlapping area of the two beam footprints,

[0182] is the area of the first beam footprint, is the area of the second beam footprint. Example:

[0183]

[0184] ​​The dashed box is the region to be detected, and the longitude and latitude of the upper left corner of the picture are about (115.625982, 38.512285), and the longitude and latitude of the lower right corner are about (125.582333, 28.672418). The adopted orbit parameter settings are shown in Table 1. Figure 4 The edge extraction result is shown in Fig. 6. Figure 5

[0185] Figure 6 The attitude result change curve diagram obtained by using matlab simulation is shown in Fig. 7. It can be seen from the result diagram that due to the protruding region in the northeast of Shandong Province, the yaw angle and the pitch angle have a trough at the 40s position. Since the roll angle is mainly to ensure that the light beam correctly irradiates the ground region, the roll angle changes of the main satellite and the auxiliary satellite are relatively consistent.

[0186] Figure 7 The beam coverage range of the attitude design result simulated by STK is shown in Fig. 8. In order to clearly reflect the long axis and the short axis of the light beam, a rectangular light beam footprint is set for simulation. Among them, the purple is the main star transmission beam, and the yellow is the auxiliary satellite transmission beam. The left figure is the light beam footprint simulation diagram based on the traditional method. The right figure is the light beam footprint simulation diagram based on the proposed method. By comparing the two figures, it can be clearly seen that the light beam footprint overlap rate of the proposed method is higher, and the effect is better.​

Claims

1. A radar attitude pointing method under a large dual-baseline angle mode, characterized by the steps of Comprise: 1) input map image, automatically identify the specified area according to color threshold; 2) input satellite orbit six numbers, UTC time and target position; 3) calculate the satellite position in WGS-84 coordinate system; 4) calculate the satellite to the target line beam vector; 5) determine the three-axis in the satellite body coordinate system; 6) calculate the attitude conversion matrix according to the principle of maximum energy, and calculate the attitude Euler angle from the attitude matrix; 7) calculate the beam footprint overlap rate.

2. The radar attitude pointing method under a large dual basis angle mode according to claim 1, characterized in that Step 1) is specifically: 1-1) read and convert image: convert from RGB color space to HSV color space; 1-2) color threshold setting: set threshold in HSV color space to identify specific area; 1-3) apply color filter: create a mask based on the set threshold to isolate the river part; 1-4) extract target area coordinates: apply the mask to the original image to extract the coordinates of the river area.

3. The radar attitude pointing method under a large dual basis angle mode according to claim 1, characterized in that The step 2) is specifically: 2-1) using six orbital elements to describe the orbital coordinate system, the six orbital elements are orbital semi-major axis ( ), orbital eccentricity ( ), orbital inclination ( ), longitude of ascending node ( ), argument of perigee ( ), true anomaly ( ). 2-2) When the satellite is in motion, the elapsed time t needs to be input to calculate the new true anomaly (θ) ).

4. The radar attitude pointing method according to claim 3, wherein Step 2-2) is specifically: 2-2-1) calculating an initial eccentricity : Known original true anomaly angle , calculate the eccentric anomaly angle : ； 2-2-2) Calculate mean anomaly : ； 2-2-3) Calculation of the orbital period : Orbital period is the time it takes for a satellite to complete one orbit around the central body, calculated by Kepler's third law: ； wherein is the standard gravitational parameter of the central celestial body; 2-2-4) Updating the Mean Anomaly : ； With time of the increase in mean anomaly, the argument of latitude will increase at a constant rate, which is determined by the orbital period . 2-2-5) Solving the Kepler equation by using iterative approximation to obtain the eccentric anomaly ; ； 2-2-6) Updating the true anomaly ; 。 5. The radar attitude pointing method under a large dual basis angle regime according to claim 1, characterized in that Step 3) is specifically: 3-1) Calculate the position of the satellite in the inertial coordinate system from the orbital six elements The calculation method is as follows: Determination of orbit plane: ； ； Calculate coordinates according to rotation matrix: ； ； ； ； 3-2) conversion from inertial coordinate system to WGS-84 coordinate system: The angle is calculated using Greenwich Mean Sidereal Time based on UTC time, and then transformed using the rotation matrix R. The conversion to WGS-84 coordinate system is calculated as follows: The rotation angle is calculated according to the Greenwich mean sidereal time calculation formula : ； Where JD is the input UTC time corresponding to the Julian day, J2000 corresponds to the Julian day 2451545, and the last two items are precession correction terms; ； 。 6. The radar attitude pointing method under a large dual basis angle regime according to claim 1, characterized in that Step 4) is specifically: 4-1) Defining the beam direction: in the satellite body coordinate system, assume the beam is transmitted along the Z axis, i.e. the beam direction is represented in the body coordinate system as ; 4-2) coordinate conversion: convert the beam direction from WGS-84 coordinate to orbit coordinate system VVLH using satellite attitude; is the inverse matrix of the coordinate axes: ； The transformation matrix around x, y, z is as follows: ； ； ； ； 。 7. The radar attitude pointing method under a large dual basis angle mode according to claim 1, characterized in that Step 5) is specifically: 5-1) Determining binary stars Axis pointing: The radar beam center always points towards the ground target particle, that is, the beam center vector of a bistatic radar in a ground-fixed system at any given time. , Represented as ； ； wherein, represents the position of satellite 1, represents the position of satellite 2, the satellite positions in the earth-fixed system can be calculated from the orbital elements and the UTC time, represents the target position; The ground-fixed system under the bistatic radar beam center vector at any time 、 Turn to the orbital coordinate system and unitize to get 、 : The inverse matrix of the coordinate axis: ； The transformation matrix around x, y, z is as follows: ； ； ， ； 5-2) Determining binary stars Axis pointing: Due to the principle of maximum energy, the dual-star beam footprints need to coincide, so the dual-star Both axes point in the direction of the velocity of the operation in the middle of the dual-star; 5-3) Determining binary stars Axis pointing: After determining axis and axis at time t can be found according to the right-hand rule axis unit vector and .

8. The radar attitude pointing method under a large dual basis angle regime according to claim 1, characterized in that Step 6) is specifically: calculate the Euler angle of the attitude generated by the converted reference coordinate system coordinates and the original body coordinate system z axis direction: wherein, referred to as a direction cosine matrix; ； For the direction cosine matrix , the ZYX Euler angles can be resolved in the following way: ； ； ； wherein here is the element of the matrix in the i-th row, j-th column.

9. The radar attitude pointing method according to claim 1, wherein Step 7) is specifically: ； wherein is the area of overlap of the two beam footprints, is the area of the first beam footprint, is the area of the second beam footprint.

Images