CN111427002B - Azimuth angle calculation method for ground measurement and control antenna pointing satellite - Google Patents

Azimuth angle calculation method for ground measurement and control antenna pointing satellite Download PDF

Info

Publication number
CN111427002B
CN111427002B CN202010197399.2A CN202010197399A CN111427002B CN 111427002 B CN111427002 B CN 111427002B CN 202010197399 A CN202010197399 A CN 202010197399A CN 111427002 B CN111427002 B CN 111427002B
Authority
CN
China
Prior art keywords
satellite
antenna
ground
station
calculation
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010197399.2A
Other languages
Chinese (zh)
Other versions
CN111427002A (en
Inventor
俞航
易灵
郑渊
刘登岭
辛阁霖
宋效正
吕旺
边志强
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shanghai Institute of Satellite Engineering
Original Assignee
Shanghai Institute of Satellite Engineering
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Shanghai Institute of Satellite Engineering filed Critical Shanghai Institute of Satellite Engineering
Priority to CN202010197399.2A priority Critical patent/CN111427002B/en
Publication of CN111427002A publication Critical patent/CN111427002A/en
Application granted granted Critical
Publication of CN111427002B publication Critical patent/CN111427002B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S1/00—Beacons or beacon systems transmitting signals having a characteristic or characteristics capable of being detected by non-directional receivers and defining directions, positions, or position lines fixed relatively to the beacon transmitters; Receivers co-operating therewith
    • G01S1/02—Beacons or beacon systems transmitting signals having a characteristic or characteristics capable of being detected by non-directional receivers and defining directions, positions, or position lines fixed relatively to the beacon transmitters; Receivers co-operating therewith using radio waves
    • G01S1/08—Systems for determining direction or position line
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Abstract

The invention discloses a method for calculating an azimuth angle of a ground measurement and control antenna pointing to a satellite, which comprises the following steps: using time of ephemeris t0、t0Number of orbital flat at time [ a ]0,e0,i0,Ω0,ω0,M0]And the groundThe geographical position information of the antenna of the area survey station is finally converted into t through the correlation calculation of the satellite orbit and the conversion calculation of a plurality of correlation coordinate systems1Pointing azimuth angle of antenna under standing center system at time: high and low angle
Figure DDA0002418106560000011
And the horizontal angle psi is used for completing the prediction process of the satellite pointing direction by the ground station antenna. The method does not depend on simulation software or excessive assumed contents, calculates the direction of the ground station antenna to the satellite by considering the actual running condition of the satellite, considers the perturbation of J2-J4 in the earth gravitational field shooting in the calculation, is suitable for long-time and high-precision satellite position prediction, can be suitable for various task requirements of the satellite in orbit, effectively solves the problem of controlling the direction of the ground station receiving antenna to the satellite, and achieves higher direction precision.

Description

Azimuth angle calculation method for ground measurement and control antenna pointing satellite
Technical Field
The invention relates to the field of satellite orbit calculation in orbit dynamics, in particular to a method for calculating an azimuth angle of a ground measurement and control antenna pointing to a satellite.
Background
The radar plays a very important role in the scientific and technological construction of China, and along with the requirements of detecting and controlling an outer space target, the practical radar is rapidly developed and is applied to a plurality of important fields such as guidance, beyond-the-horizon detection and the like at present. With the development of the space detection technology, higher and higher requirements are put forward on the tracking and searching capabilities of a radar antenna, because a satellite signal is weak and has strong directivity, in order to capture a communication signal on a moving satellite, the deviation between the attitude of the antenna and the position of the satellite must be adjusted in real time to meet the communication requirement, because the signal-to-noise ratio of link transmission information is reduced due to the satellite-to-ground directional deviation, and if the signal-to-noise ratio exceeds the maximum station tolerance, the signal loss phenomenon can even occur. This requires that the radar antenna must adjust the pointing direction according to the command to track the moving target in real time. Therefore, the dynamic precision of the radar antenna pointing process becomes one of the important indexes of the antenna system function, and the design of the pointing calculation method with high pointing precision has general practical significance.
The method for forecasting the pointing direction of the ground survey station antenna to the satellite mainly comprises the steps of calculating an antenna pointing azimuth angle at a target moment according to orbit information and time information of the satellite and position information of the ground survey station antenna, and realizing accurate pointing to the satellite at the target moment through antenna pointing control. The prediction of the position of a satellite at a given moment has become an increasingly important issue in recent years. The high-precision orbit prediction is an important technology in the aerospace technology, plays an important role in satellite orbit design and orbit optimization, and can provide reliable orbit information reference for satellite in-orbit tasks such as antenna pointing, tracking and positioning and the like.
The existing satellite-ground pointing algorithm research in China mostly focuses on the optimization design of the pointing of a satellite to a ground survey station under the condition that the position of the ground survey station is fixed, and the research on the pointing of a ground survey station antenna to the satellite is less. Under the condition that the ground survey station is movable, under the condition that the geographical longitude, latitude and elevation of the ground survey station are known and the track information of the on-orbit aircraft is tracked, the pointing control of the aircraft is automatically completed, and the pointing azimuth angle of the antenna is calculated in real time. At present, a satellite orbit recursion method and a high-precision orbit determination algorithm applied to a ground station are generally realized by a high-performance computer, wherein the high-precision computer comprises a high-precision integral algorithm and a high-precision dynamic model. Aiming at the actual situation, the invention provides a ground measurement station real-time positioning method with higher precision for a ground measurement station antenna, and the method can realize the prediction of the autonomous pointing direction of the satellite by the ground measurement station antenna.
The patent "design method of deep space probe antenna pointing" (patent number: CN104369877A) describes a method of pointing a deep space probe antenna to the ground center, which is used to realize the deep space probe antenna to the ground center orientation. The patent is directed to the orientation of the antenna to the geocenter and not to the given position of the earth surface, and the pointing vector of the detector antenna to the geocenter is directly given, and no algorithm for calculating the satellite position through the orbit parameters exists. The method is different from the method in that a method for calculating the satellite pointing direction of the ground station aiming at the given position of the earth surface is designed, the positioning calculation of the ground station is completed, and a calculation process for calculating the ground station-satellite pointing vector through the satellite orbit parameters at the given moment is designed.
The patent "a simulation analysis method of pointing angle of data transmission antenna" (patent No. CN105184002A) introduces a method for calculating pointing direction of satellite-borne data transmission antenna to ground station, which uses the existing satellite orbit simulation software STK to perform simulation solution on the actual position of the satellite and calculates the two-dimensional pointing angle of the data transmission antenna. The disadvantage of the patent is that the satellite position calculation depends on the satellite orbit simulation software STK, no specific calculation process is needed, the description of the coordinate system conversion is simple, and no algorithm of a conversion matrix is given. The invention has the advantages of providing a method for calculating the actual position of the satellite according to the orbit parameters of the satellite at the appointed time without depending on STK software and designing a set of detailed calculation flow of a related coordinate system conversion matrix.
The patent "a method for controlling the pointing direction of a dual-axis antenna to the ground around a moon satellite" (patent number: CN101204994A) describes a method for calculating the pointing direction of a satellite to the earth center around a moon satellite, which estimates the position of the satellite according to ephemeris data on the ground, calculates the visible area of the satellite to the earth, and calculates the pointing angle of the dual-axis antenna. The patent is directed to the geocenter, does not orient the surface location, and is mainly calculated in combination with a moon-related coordinate system. The invention is different from the method in that the calculation is mainly combined with the earth and the earth surface position related coordinate system to complete the position calculation of the ground survey station antenna and the directional calculation of the ground survey station antenna to the satellite, and the definition and the calculation method of the antenna directional angle are different.
Li Dan et al in the article "Low orbit satellite orbit prediction algorithm based on orbit number" (see optical precision engineering, 2016, 10) propose a method for predicting satellite orbits by using elliptic curves, but in the solution process, partial differentiation of coefficients needs to be calculated.
In the Chinese invention patent 'an on-satellite autonomous orbit extrapolation method suitable for a circular orbit satellite' (patent number: CN103995800A), a method for orbit recursion suitable for the circular orbit satellite is introduced. But this method only takes perturbation into account.
Based on the above consideration, the azimuth angle calculation method for the ground measurement and control antenna pointing satellite disclosed by the invention has the advantages of high precision and suitability for long-term prediction.
Disclosure of Invention
Aiming at the defects in the prior art, the invention aims to provide a method for calculating the azimuth angle of a ground measurement and control antenna pointing to a satellite.
The method for calculating the azimuth angle of the ground measurement and control antenna pointing to the satellite comprises the following steps:
using ephemeris time t0、t0Number of orbital flat at time [ a ]0,e0,i0,Ω0,ω0,M0]And the geographical position information of the ground survey station antenna is finally converted into t through the related calculation of the satellite orbit and the conversion calculation of a plurality of related coordinate systems1The antenna pointing azimuth angle at the moment under the standing heart system is as follows: high and low angle
Figure BDA0002418106540000031
And the horizontal angle psi is used for completing the prediction process of the satellite pointing direction by the ground station antenna. Further, ephemeris time t0、t0Number of orbital plateau at time [ a ]0,e0,i0,Ω0,ω0,M0]And the geographic position information of the ground station antenna are known quantities and are directly obtained by inquiring the ephemeris of the satellite.
Preferably, the method comprises the following steps:
the input parameter non-singular processing steps:
when the track eccentricity is small, i.e. near circular track, e ≈ 0, to avoid the occurrence of singularities in the calculation, we turn to 3 singularity-free variables, i.e.:
ξ0＝e0cos(ω0)
η0＝-e0sin(ω0)
λ0＝ω0+M0
preferably, the method further comprises the following steps:
input parameter normalization processing:
normalization processing is performed for the recursion time dt:
dt＝t1-t0
Figure BDA0002418106540000032
wherein, the first and the second end of the pipe are connected with each other,
t0recursion of the starting moment for the track;
t1representing the track recursion end time;
dtnthe result of normalization processing of recursion time is shown, and subscript n shows normalization;
for semi-major axis a0And (3) carrying out normalization treatment:
Figure BDA0002418106540000033
the normalized unit of the time is
Figure BDA0002418106540000041
Wherein the content of the first and second substances,
ge is a gravitational constant;
re represents the earth equatorial radius.
Preferably, the method further comprises the following steps:
perturbation item calculation step:
considering the perturbations of items J2 to J4 in the earth gravity field, including the calculation of first-order long-term, first-order short-term and second-order long-term terms, the expressions for the individual perturbation terms are as follows:
first order long term calculation:
Figure BDA0002418106540000042
Figure BDA0002418106540000043
Figure BDA0002418106540000044
wherein the content of the first and second substances,
Ω1representing the first-order long period variable quantity of the right ascension of the satellite orbit intersection point;
ω1representing the first-order long period variation of the argument of the satellite orbit in the near place;
λ1represents the first-order long-period variation of the singularity-free variable λ introduced by the calculations herein;
p is a radius, and the expression is p ═ a × (1-e)0 2) Mean angular velocity of track
Figure BDA0002418106540000045
J2＝1.624×10-3；
First order short period term calculation:
Figure BDA0002418106540000046
Figure BDA0002418106540000047
Figure BDA0002418106540000048
Figure BDA0002418106540000049
Figure BDA00024181065400000410
Figure BDA0002418106540000051
wherein the content of the first and second substances,
Figure BDA0002418106540000052
representing the first-order short period variable quantity of the semi-major axis of the satellite orbit;
Figure BDA0002418106540000053
representing the first-order short-period variation of the satellite orbit inclination angle;
Figure BDA0002418106540000054
the first-order short-period variation of the right ascension of the satellite orbit intersection point is represented;
Figure BDA0002418106540000055
respectively representing the first-order short-period variable quantities of three singularity-free variables xi, eta and lambda introduced by the calculation of the text;
the u is calculated by a first-order long-term, and the calculation process is as follows:
ξz1＝ξ0 cos(ω1 dtn)+η0 sin(ω1 dtn)
ηz1＝η0 cos(ω1 dtn)-ξ0 sin(ω1 dtn)
λz1＝λ0+(n+λ1)dtn
Figure BDA0002418106540000056
ωz1＝arc cos(ξz1/ez1)
Mz1＝mod(λz1-ωz1，2π)
Figure BDA0002418106540000057
u＝fz1+ωz1
wherein the content of the first and second substances,
ξz1、ηz1、λz1respectively representing the integration results of three singularity-free variables xi, eta and lambda introduced by the calculation of the text according to the first-order long-term variation;
ez1an integral result which represents the eccentricity of the satellite orbit and takes the first-order long period variable quantity as an integral quantity;
ωz1an integral result which represents the argument of the satellite orbit near place and takes the first-order long period variable quantity as integral quantity;
Mz1an integration result which represents the mean-near point angle of the satellite orbit and takes the first-order long period variable quantity as an integral quantity;
fz1an integral result which represents the true near point angle of the satellite orbit and takes the first-order long period variable quantity as an integral quantity;
second-order long-term calculation:
Figure BDA0002418106540000061
Figure BDA0002418106540000062
Figure BDA0002418106540000063
wherein the content of the first and second substances,
Ω2representing the second-order long-period variation of the right ascension of the satellite orbit;
ξ2、λ2respectively representing second-order long-period variable quantities of singularity-free variables xi and lambda introduced by calculation in the text;
J3＝2.5356×10-6,J4＝7.1022×10-6。
preferably, the method further comprises the following steps:
and (3) calculating a recursion main formula:
for introduced singularity-free variables, an analytic solution structure is carried out by combining a long-term item perturbation item and a short-period item perturbation item, the long-period item perturbation is ignored, and a track recursion main formula is as follows:
Figure BDA0002418106540000064
Figure BDA0002418106540000065
Figure BDA0002418106540000066
Figure BDA0002418106540000067
Figure BDA0002418106540000068
Figure BDA0002418106540000069
wherein the content of the first and second substances,
αtrepresents t1At the moment, the normalization result of the satellite orbit semi-major axis;
isrepresents t1At that time, the satellite orbit inclination;
Ωsrepresents t1At the moment, the rising point of the satellite orbit is right ascension;
ξs、ηs、λsrespectively represent t1At the moment, the values of the three introduced singularity-free variables ξ, η, λ are calculated.
Preferably, the method further comprises the following steps:
a singular point-free variable reduction step:
after the calculation is finished, 3 singularity-free variables are restored
Figure BDA0002418106540000071
ωs＝arc cos(ξs/es)
Ms＝mod(λs-ωs,2π)
Preferably, the method further comprises the following steps:
and (3) normalization variable reduction step:
will be the long axis a of the satellite orbittReduction to conventional units:
as＝at×Re
a is describedsUnit: and m is selected.
Preferably, the method further comprises the following steps:
calculating the position of the inertial system satellite:
input t1Instantaneous number of satellite orbit [ a ] at times,es,is,Ωs,ωs,Ms]Outputting the component R of the position vector of the satellite in the J2000.0 coordinate systemwECIThe calculation method comprises the following steps:
RwECI＝Q*rp
wherein the rotation matrix Q is described in a 3-1-3 rotation order:
Figure BDA0002418106540000072
vector rp：
Figure BDA0002418106540000073
Wherein M is1True proximal angle:
Figure BDA0002418106540000074
preferably, the method further comprises the following steps:
and a step of calculating the position of the earth-fixed satellite:
inputting a target time t1And said calculated position R of the satellite in the inertial systemwECIOutput t of1Position R of time satellite under earth's fixationwECFThe specific method comprises the following steps:
according to a given target time t1Calculating a second counting value t from epoch J2000.0 (1/12/2000) to a predetermined target timecInputting t1Year, month, day, hour, minute, second of the moment, calculate julian day JD:
Figure BDA0002418106540000081
wherein, floor () is a round-down operation;
calculating a second counting value t from epoch J2000.0 to a given target time according to the Confucian day JDc：
tc＝(JD-2455197.5)×86400+315547200
According to the calculated epoch J2000.0 to the second counting value t of the given target timecCalculating a rotation matrix ER, a nutation matrix NR and a precision matrix PR of the earth, and calculating a conversion matrix M from an inertia system to a ground-fixed systemECI2ECF：
MECI2ECF＝ER*NR*PR
And according to t obtained by the calculation1Position R of time satellite under inertial systemwECICalculating t1Position R of time satellite under earth's fixationwECF：
RwECF＝MECI2ECF*RwECI
Preferably, the method further comprises the following steps:
and the antenna position calculation step of the earth-fixed system survey station:
according to the longitude, latitude and elevation of the given ground measurement station antenna, the position R of the ground measurement station antenna under the ground fixation system is calculatedtECFThe specific method comprises the following steps:
inputting longitude lon, latitude lat and elevation h of an antenna of the ground station;
computing coordinate components G1, G2:
Figure BDA0002418106540000082
Figure BDA0002418106540000083
wherein f is the geometric oblateness of an earth ellipsoid, and f is 1/298.257;
calculating the position R of the ground survey station antenna under the ground fixing systemtECF：
Figure BDA0002418106540000091
Preferably, the method further comprises the following steps:
and a step of calculating the position of the satellite in the station center system:
inputting calculated t1Position R of time satellite under earth's fixationwECFAnd said calculated position R of the ground station antenna under the ground anchor systemtECFOutputting the position R of the satellite under the station center systemwCTThe specific method comprises the following steps:
under the system of the station center, a conversion matrix M from the earth fixation system to the system of the station center is calculatedECF2CTDescribed as one rotation about the Z-axis of the earth fixation system and one rotation about the X-axis of the earth fixation system:
MECF2CT＝Rx(90°-lat)Rz(90°+lon)
wherein lon is the geographic longitude of the antenna of the ground survey station; lat is the geographical latitude of the antenna of the ground station;
Figure BDA0002418106540000092
Figure BDA0002418106540000093
and according to said t1Position R of time satellite under earth's fixationwECFPosition R of the ground station antenna under the ground anchortECFTranslating the origin of coordinates from the geocentric to the antenna of the ground survey station, and calculating t1Position R of time satellite under the system of the center of the stationwCT：
RwCT＝MECF2CT*(RwECF-RtECF)
Preferably, the method further comprises the following steps:
the method comprises the following steps of:
inputting the position of the satellite under the station center system, and outputting the pointing azimuth angle of the ground station antenna at the time t 1: high low angle
Figure BDA0002418106540000094
The horizontal angle psi is specifically as follows:
the high and low angles and the horizontal angle are defined under the station center system OCTXCTYCTZCTIn, OCTRepresenting the origin of the coordinate system, XCTX-axis, Y, representing a coordinate systemCTY-axis, O, representing a coordinate systemCTXCTYCTZCTThe plane, elevation angle, of X-axis and Y-axis of the coordinate system
Figure BDA0002418106540000101
To point to a vector RwCTAnd OCTXCTYCTAngle of plane, define RwCTVector sum of OCTZCTThe included angle is less than 90 degrees and is positive; the horizontal angle psi being the director vector RwCTAt OCTXCTYCTProjection of plane and OCTXCTAngle of axis, defined around OCTZCTShaft driven OCTXCTShaft clockwise steering pointing vector RwCTAt OCTXCTYCTThe projection of the surface is positive, and the antenna pointing azimuth angle is found according to the definition. Assuming that the position of the ground survey station is located at the origin of the station center system, the projection of the antenna pointing vector of the survey station under the station center system is RwCTRecord of RwCTComprises the following steps:
Figure BDA0002418106540000102
wherein, the first and the second end of the pipe are connected with each other,
xCT、yCT、zCTrespectively representing the pointing vectors R of antennas of the stationswCTAt the station center is OCTXCTYCTZCTCoordinate components of corresponding X-axis, Y-axis and Z-axis;
calculation of elevation angle of antenna pointing azimuth
Figure BDA0002418106540000103
The horizontal angle ψ is:
Figure BDA0002418106540000104
Figure BDA0002418106540000105
and according to xCT、yCT、zCTPositive and negative of (2), angle of elevation
Figure BDA0002418106540000106
Dividing the horizontal angle psi into corresponding angle ranges to finish the process of forecasting the antenna pointing direction:
if z isCTMore than or equal to 0, the height and angle will be changed
Figure BDA0002418106540000107
Division into ranges
Figure BDA0002418106540000108
If xCT≥0,yCTThe horizontal angle psi is divided into ranges
Figure BDA0002418106540000109
If xCT＜0,yCTThe horizontal angle psi is divided into ranges
Figure BDA00024181065400001010
If xCT＜0,yCT< 0, divide the horizontal angle psi into the range
Figure BDA00024181065400001011
If xCT≥0,yCTIf < 0, the horizontal angle psi is divided into ranges
Figure BDA00024181065400001012
Compared with the prior art, the invention has the following beneficial effects:
the method does not depend on simulation software or excessive assumed contents, considers the actual operation condition of the satellite to calculate the direction of the ground station antenna to the satellite, considers the perturbation from J2 to J4 in the earth gravity field shooting in the calculation, is suitable for long-time and high-precision satellite position prediction, can be suitable for various task requirements of the satellite in orbit, effectively solves the problem of controlling the direction of the ground station receiving antenna to the satellite, and achieves higher direction precision.
Drawings
Other features, objects and advantages of the invention will become more apparent upon reading of the detailed description of non-limiting embodiments with reference to the following drawings:
fig. 1 is a schematic flow chart of a method for calculating an azimuth angle of a ground measurement and control antenna pointing to a satellite.
Fig. 2 is a schematic diagram of autonomous prediction of the satellite orbit position.
FIG. 3 is a schematic diagram of the pointing direction of the ground station antenna to the satellite.
FIG. 4 is a view of the station center system OCTXCTYCTZCTSchematic representation.
FIG. 5 shows the elevation angle of the standing heart system
Figure BDA0002418106540000112
And horizontal angle psi.
Detailed Description
The present invention will be described in detail with reference to specific examples. The following examples will assist those skilled in the art in further understanding the invention, but are not intended to limit the invention in any way. It should be noted that it would be obvious to those skilled in the art that various changes and modifications can be made without departing from the spirit of the invention. All falling within the scope of the present invention.
The method for calculating the azimuth angle of the ground measurement and control antenna pointing to the satellite comprises the following steps:
using ephemeris time t0、t0Number of orbital flat at time [ a ]0,e0,i0,Ω0,ω0,M0]And the geographic position information of the ground survey station antenna through the correlation calculation of the satellite orbitAnd converting calculation of multiple related coordinate systems into t1Pointing azimuth angle of antenna under standing center system at time: high and low angle
Figure BDA0002418106540000111
And the horizontal angle psi is used for completing the prediction process of the satellite pointing direction by the ground station antenna. Further, ephemeris time t0、t0Number of orbital plateau at time [ a ]0,e0,i0,Ω0,ω0,M0]And the geographic position information of the ground station antenna are known quantities and are directly obtained by inquiring the ephemeris of the satellite.
Specifically, the method comprises the following steps:
the input parameter non-singular processing steps:
when the track eccentricity is small, i.e. near circular track, e ≈ 0, to avoid the occurrence of singularities in the calculation, we turn to 3 singularity-free variables, i.e.:
ξ0＝e0 cos(ω0)
η0＝-e0 sin(ω0)
λ0＝ω0+M0
specifically, the method further comprises the following steps:
input parameter normalization processing step:
normalization processing is performed for the recursion time dt:
dt＝t1-t0
Figure BDA0002418106540000121
wherein the content of the first and second substances,
t0recursion of the starting time for the track;
t1representing the track recursion end time;
dtnthe result of normalization processing of the recursion time is shown, and the subscript n shows normalization;
for semi-major axis a0Go on to returnNormalization treatment:
Figure BDA0002418106540000122
the normalized unit of the time is
Figure BDA0002418106540000123
Wherein the content of the first and second substances,
ge is a gravitational constant;
re represents the earth equatorial radius.
Specifically, the method further comprises the following steps:
perturbation item calculation step:
considering the perturbations of items J2 to J4 in the earth gravity field, including the calculation of first-order long-term, first-order short-term and second-order long-term terms, the expressions for the individual perturbation terms are as follows:
first order long term calculation:
Figure BDA0002418106540000131
Figure BDA0002418106540000132
Figure BDA0002418106540000133
wherein the content of the first and second substances,
Ω1representing the first-order long period variable quantity of the right ascension of the satellite orbit intersection point;
ω1representing the first-order long-period variation of the argument of the satellite orbit in the near place;
λ1represents the first-order long-period variation of the singularity-free variable λ introduced by the calculations herein;
p is a half-diameter and the expression isp＝a×(1-e0 2) Mean angular velocity of track
Figure BDA0002418106540000134
J2＝1.624×10-3；
First order short period term calculation:
Figure BDA0002418106540000135
Figure BDA0002418106540000136
Figure BDA0002418106540000137
Figure BDA0002418106540000138
Figure BDA0002418106540000139
Figure BDA00024181065400001310
wherein the content of the first and second substances,
Figure BDA00024181065400001311
representing the first-order short period variable quantity of the semi-major axis of the satellite orbit;
Figure BDA00024181065400001312
representing the first-order short-period variation of the satellite orbit inclination angle;
Figure BDA00024181065400001313
the first-order short-period variation of the right ascension of the satellite orbit intersection point is represented;
Figure BDA00024181065400001314
respectively representing first-order short-period variable quantities of three singularity-free variables xi, eta and lambda introduced by the calculation of the text;
the u is calculated by a first-order long-term, and the calculation process is as follows:
ξz1＝ξ0 cos(ω1 dtn)+η0 sin(ω1 dtn)
ηz1＝η0 cos(ω1 dtn)-ξ0 sin(ω1 dtn)
λz1＝λ0+(n+λ1)dtn
Figure BDA0002418106540000141
ωz1＝arc cos(ξz1/ez1)
Mz1＝mod(λz1-ωz1，2π)
Figure BDA0002418106540000142
u＝fz1+ωz1
wherein the content of the first and second substances,
ξz1、ηz1、λz1respectively representing the integral results of three singularity-free variables xi, eta and lambda introduced by the calculation of the text according to the first-order long-term variation quantity;
ez1an integral result which represents the eccentricity of the satellite orbit and takes the first-order long period variable quantity as an integral quantity;
ωz1representing variation of amplitude angle of satellite orbit in first-order long periodThe chemical quantity is an integral result of the integral quantity;
Mz1an integration result which represents the mean-near point angle of the satellite orbit and takes the first-order long period variable quantity as an integral quantity;
fz1an integral result which represents the true near point angle of the satellite orbit and takes the first-order long period variable quantity as integral quantity;
second-order long-term calculation:
Figure BDA0002418106540000143
Figure BDA0002418106540000144
Figure BDA0002418106540000145
Figure BDA0002418106540000146
wherein the content of the first and second substances,
Ω2representing the second-order long-period variation of the right ascension of the satellite orbit;
ξ2、λ2respectively representing second-order long-period variable quantities of singularity-free variables xi and lambda introduced by calculation;
J3＝2.5356×10-6,J4＝7.1022×10-6。
specifically, the method further comprises the following steps:
and (3) calculating a recursion main formula:
for introduced singularity-free variables, an analytic solution structure is carried out by combining a long-term item perturbation item and a short-period item perturbation item, the long-period item perturbation is ignored, and a track recursion main formula is as follows:
Figure BDA0002418106540000151
Figure BDA0002418106540000152
Figure BDA0002418106540000153
Figure BDA0002418106540000154
Figure BDA0002418106540000155
Figure BDA0002418106540000156
wherein the content of the first and second substances,
αtrepresents t1At the moment, the normalization result of the satellite orbit semi-major axis;
isrepresents t1At that time, the satellite orbit inclination;
Ωsdenotes t1At the moment, the rising point of the satellite orbit is right ascension;
ξs、ηs、λsrespectively represent t1At the moment, the values of the three introduced singularity-free variables ξ, η, λ are calculated.
Specifically, the method further comprises the following steps:
a singular point-free variable reduction step:
after the calculation is finished, 3 singularity-free variables are restored
Figure BDA0002418106540000157
ωs＝arc cos(ξs/es)
Ms＝mod(λs-ωs,2π)
Specifically, the method further comprises the following steps:
and (3) normalization variable reduction step:
will be the long axis a of the satellite orbittReduction to conventional units:
as＝at×Re
a is as describedsUnit: and m is selected.
Specifically, the method further comprises the following steps:
calculating the position of the inertial system satellite:
input t1Instantaneous number of satellite orbits [ a ] of times,es,is,Ωs,ωs,Ms]Outputting the component R of the position vector of the satellite in the J2000.0 coordinate systemwECIThe calculation method comprises the following steps:
RwECI＝Q*rp
wherein the rotation matrix Q is described in a 3-1-3 rotation order:
Figure BDA0002418106540000161
vector rp：
Figure BDA0002418106540000162
Wherein, M1True proximal angle:
Figure BDA0002418106540000163
specifically, the method further comprises the following steps:
and a step of calculating the position of the earth-fixed satellite:
inputting a target time t1And said calculated position R of the satellite in the inertial systemwECIOutput t of1Position R of time satellite under earth's fixationwECFThe specific method comprises the following steps:
according to a given target time t1Calculating a second counting value t from epoch J2000.0 (1/12/2000) to a predetermined target timecInputting t1Year, month, day, hour, minute, second of the moment, calculate julian day JD:
Figure BDA0002418106540000164
wherein, floor () is a rounding-down operation;
calculating a second count value t from epoch J2000.0 to a given target time based on the julian day JDc：
tc＝(JD-2455197.5)×86400+315547200
According to the second counting value t from the epoch J2000.0 to the given target timecCalculating an earth rotation matrix ER, a nutation matrix NR and a time difference matrix PR, and calculating a conversion matrix M from an inertia system to a ground-fixed systemECI2ECF：
MECI2ECF＝ER*NR*PR
And according to t obtained by the calculation1Position R of time satellite under inertial systemwECICalculating t1Position R of time satellite under earth's fixationwECF：
RwECF＝MECI2ECF*RwECI
Specifically, the method further comprises the following steps:
and the antenna position calculation step of the earth-fixed system survey station:
according to the longitude, latitude and elevation of the given ground measurement station antenna, the position R of the ground measurement station antenna under the ground fixation system is calculatedtECFThe specific method comprises the following steps:
inputting longitude lon, latitude lat and elevation h of an antenna of the ground station;
computing coordinate components G1, G2:
Figure BDA0002418106540000171
Figure BDA0002418106540000172
wherein f is the geometric oblateness of an earth ellipsoid, and f is 1/298.257;
calculating the position R of the ground survey station antenna under the ground fixing systemtECF：
Figure BDA0002418106540000173
Specifically, the method further comprises the following steps:
and a step of calculating the position of the station center system satellite:
inputting calculated t1Position R of time satellite under earth fixed systemwECFAnd said calculated position R of the ground station antenna under the ground anchor systemtECFOutputting the position R of the satellite under the station center systemwCTThe specific method comprises the following steps:
under the system of the station center, a conversion matrix M from the earth fixation system to the system of the station center is calculatedECF2CTDescribed as one rotation about the Z axis of the earth fixation and one rotation about the X axis of the earth fixation:
MECF2CT＝Rx(90°-lat)Rz(90°+lon)
wherein lon is the geographic longitude of the antenna of the ground survey station; lat is the geographical latitude of the antenna of the ground station;
Figure BDA0002418106540000181
Figure BDA0002418106540000182
and according to said t1Position R of time satellite under earth fixed systemwECFPosition R of the ground station antenna under the ground anchortECFTranslating the origin of coordinates from the geocenter to the antenna of the ground survey station, and calculating t1Time of day satellite in-station systemPosition R ofwCT：
RwCT＝MECF2CT*(RwECF-RtECF)
Specifically, the method further comprises the following steps:
the method comprises the following steps of:
inputting the position of the satellite under the station center system and outputting t1The antenna pointing azimuth angle of the ground survey station at any moment: high low angle
Figure BDA0002418106540000183
The horizontal angle psi is specifically as follows:
the high and low angles and the horizontal angle are defined under the station center system OCTXCTYCTZCTIn, OCTRepresenting the origin of the coordinate system, XCTX-axis, Y, representing a coordinate systemCTY-axis, O, representing a coordinate systemCTXCTYCTThe plane, elevation angle, of X-axis and Y-axis of the coordinate system
Figure BDA0002418106540000184
To point to a vector RwCTAnd OCTXCTYCTAngle of plane, define RwCTVector sum of OCTZCTThe included angle is less than 90 degrees and is positive; the horizontal angle psi being the director vector RwCTAt OCTXCTYCTProjection of plane and OCTXCTAngle of axis, defined around OCTZCTShaft driven OCTXCTShaft clockwise steering pointing vector RwCTAt OCTXCTYCTThe projection of the surface is positive, and the antenna pointing azimuth angle is found according to the definition. Assuming that the position of the ground survey station is located at the origin of the station center system, the projection of the antenna pointing vector of the survey station under the station center system is RwCTRecord RwCTComprises the following steps:
Figure BDA0002418106540000185
wherein, the first and the second end of the pipe are connected with each other,
xCT、yCT、zCTrespectively representing the pointing vectors R of antennas of the stationswCTAt the station center system OCTXCTYCTZCTCoordinate components of corresponding X-axis, Y-axis and Z-axis;
calculation of elevation angle of antenna pointing azimuth
Figure BDA0002418106540000191
The horizontal angle ψ is:
Figure BDA0002418106540000192
Figure BDA0002418106540000193
and according to xCT、yCT、zCTPositive and negative of (2), angle of elevation
Figure BDA0002418106540000194
Dividing the horizontal angle psi into corresponding angle ranges to finish the process of forecasting the antenna pointing direction:
if z isCTMore than or equal to 0, the height and angle will be changed
Figure BDA0002418106540000195
Division into ranges
Figure BDA0002418106540000196
If xCT≥0,yCTThe horizontal angle psi is divided into ranges
Figure BDA0002418106540000197
If xCT＜0,yCTThe horizontal angle psi is divided into ranges
Figure BDA0002418106540000198
If xCT＜0,yCTIf < 0, the horizontal angle psi is divided into ranges
Figure BDA0002418106540000199
If xCT≥0,yCT< 0, divide the horizontal angle psi into the range
Figure BDA00024181065400001910
The present invention will be described more specifically below with reference to preferred examples.
Preferred example 1:
the technical problem to be solved by the invention is as follows: the method comprises the steps of performing satellite orbit correlation calculation and conversion calculation of a plurality of correlation coordinate systems by satellite ephemeris data and geographical position information of a ground measurement station antenna at a given target moment, and finally converting the satellite ephemeris data and the geographical position information into an antenna pointing azimuth angle under a station center system to finish the prediction process of the ground measurement station antenna on the satellite azimuth.
The invention combines a practical engineering situation: under the condition that the ground survey station is movable, under the condition that the geographical latitude, longitude and altitude of the ground survey station are known and the track information of the on-orbit aircraft is tracked, the pointing tracking of the aircraft is automatically completed, and the pointing azimuth angle of the antenna is calculated in real time. The ground measurement station real-time positioning method is used for the ground measurement station antenna with high precision, and the autonomous pointing direction of the satellite can be forecasted by the ground measurement station antenna.
The method for calculating the azimuth angle of the ground measurement and control antenna pointing to the satellite calculates the position of the satellite in real time through satellite ephemeris information, considers influence factors of a coordinate system conversion relation comprehensively, has high calculation precision, provides the pointing azimuth angle definition suitable for the ground measurement station antenna, and effectively meets the prediction requirement of the ground measurement station antenna on the satellite pointing azimuth in real time.
In order to realize the method, the following technical scheme is adopted:
a method for calculating the azimuth angle of a ground measurement and control antenna pointing to a satellite comprises the following steps:
1. input parameter non-singular processing module
When the track eccentricity is very small (near-circular track, e ≈ 0), in order to avoid the occurrence of singular points in the calculation, the method is converted into 3 singular point-free variables, namely
ξ0＝e0 cos(ω0)
η0＝-e0 sin(ω0)
λ0＝ω0+M0
2. Input parameter normalization processing module
The normalization process is performed with respect to the recursion time dt,
dt＝t1-t0
Figure BDA0002418106540000201
for semi-major axis a0The normalization treatment is carried out, and the normalization treatment is carried out,
Figure BDA0002418106540000202
the normalized unit of the time is that,
Figure BDA0002418106540000203
wherein Ge is a gravitational constant;
the normalized unit of the length is Re-6378140 m
3. Perturbation item calculation module
Considering the perturbations of items J2 to J4 in the earth gravity field, including computing first order long term, first order short term, and second order long term terms, the expressions for the respective perturbation terms:
first order long term calculation:
Figure BDA0002418106540000211
Figure BDA0002418106540000212
Figure BDA0002418106540000213
the p is a radius, and the expression is p ═ a × (1-e)0 2) Mean angular velocity of track
Figure BDA0002418106540000214
J2＝1.624×10-3。
First order short period term calculation:
Figure BDA00024181065400002113
Figure BDA0002418106540000215
Figure BDA0002418106540000216
Figure BDA0002418106540000217
Figure BDA0002418106540000218
Figure BDA0002418106540000219
and u is calculated by a first-order long-term.
Figure BDA00024181065400002110
Figure BDA00024181065400002111
ξz1＝ξ0 cos(ωc1dtn)+η0 sin(ωc1dtn)
ηz1＝η0 cos(ωc1dtn)-ξ0 sin(ωc1dtn)
λz1＝λ0+(n+λc1)dtn
Figure BDA00024181065400002112
ωz1＝arc cos(ξz1/ez1)
Mz1＝mod(λz1-ωz1,2π)
Figure BDA0002418106540000221
u＝fz1+ωz1
Second-order long-term calculation:
Figure BDA0002418106540000222
Figure BDA0002418106540000223
Figure BDA0002418106540000224
said, J3＝2.5356×10-6,J4＝7.1022×10-6。
4. Recursion main formula calculation module
And for the introduced singularity-free variable, carrying out analytical solution construction by combining a long-term item perturbation item and a short-period item perturbation item, and ignoring the long-period item perturbation. The main formula of the track recursion is as follows:
Figure BDA0002418106540000225
Figure BDA0002418106540000226
Figure BDA0002418106540000227
Figure BDA0002418106540000228
Figure BDA0002418106540000229
Figure BDA00024181065400002210
5. variable recovery module without singularity
After the calculation is finished, 3 singularity-free variables are restored
Figure BDA00024181065400002211
ωs＝arc cos(ξs/es)
Ms＝mod(λs-ωs,2π)
6. Normalization variable reduction module
Will the satellite orbit major axis atReduction to conventional units:
as＝at×Re
a is describedsUnit: and m is selected.
7. Inertial system satellite position module
Passing through t1Instantaneous number of satellite orbit [ a ] at times,es,is,Ωs,ωs,Ms]Calculating the component R of the position vector of the satellite in the J2000.0 coordinate systemwECIThe calculation method comprises the following steps:
RwECI＝Q*rp
wherein the rotation matrix Q is described in a 3-1-3 rotation order:
Figure BDA0002418106540000231
vector rp：
Figure BDA0002418106540000232
Wherein M is1True proximal angle:
Figure BDA0002418106540000233
8. earth fixed satellite position module
According to a given target time t1And the position R of the satellite calculated in the step (7) in the inertial systemwECICalculating t1Position R of time satellite under earth's fixationwECFThe method comprises the following steps:
according to a given target time t1Calculating a second counting value t from epoch J2000.0 (1/12/2000) to a predetermined target timecInputting t1Time of dayYear (year), month (month), day (day), hour (hour), minute (min), second (sec) of (UTC time), julian day JD is calculated:
Figure BDA0002418106540000234
wherein floor () is a round-down operation.
Calculating a second counting value t from epoch J2000.0 (1 month, 1 day, 12 hours in 2000) to a given target time according to the julian day JDc：
tc＝(JD-2455197.5)×86400+315547200
A second counting value t from the epoch J2000.0 (1 month, 1 day, 12 hours in 2000) to a given target timecAnd calculating a terrestrial rotation matrix ER, a nutation matrix NR and a time offset matrix PR, wherein the term is not considered in the invention because polar shift has little influence on the calculation of the conversion matrix. Calculating a transformation matrix M from an inertial system to a ground-based systemECI2ECF：
MECI2ECF＝ER*NR*PR
And calculating t according to the step (7)1Position R of time satellite under inertial systemwECICalculating t1Position R of time satellite under earth's fixationwECF：
RwECF＝MECI2ECF*RwECI
9. Antenna position module of ground fixing system measuring station
The method is characterized in that the position R of the antenna of the ground measuring station under the ground fixation system is calculated according to the longitude, latitude and elevation of the given antenna of the ground measuring stationtECFThe method comprises the following steps:
and inputting longitude lon, latitude lat and elevation h of the antenna of the ground station.
Computing coordinate components G1, G2:
Figure BDA0002418106540000241
Figure BDA0002418106540000242
wherein f is the geometric oblateness of the earth ellipsoid, and f is 1/298.257.
Calculating the position R of the ground survey station antenna under the ground fixing systemtECF：
Figure BDA0002418106540000243
10. Satellite position module for station center system
T calculated according to the step (8)1Position R of time satellite under earth's fixationwECFAnd the position R of the ground survey station antenna under the ground fixed system calculated in the step (9)tECFCalculating t1Position R of time satellite under the system of the center of the stationwCTThe method comprises the following steps:
under the system of the station center, a conversion matrix M from the earth fixation system to the system of the station center is calculatedECF2CTDescribed as one rotation about the Z-axis of the earth fixation system and one rotation about the X-axis of the earth fixation system:
MECF2CT＝Rx(90°-lat)Rz(90°+lon)
wherein lon is the geographic longitude of the antenna of the ground survey station; lat is the geographic latitude of the ground station antenna.
Figure BDA0002418106540000251
Figure BDA0002418106540000252
And according to t calculated in the step (8)1Position R of time satellite under earth's fixationwECFThe position R of the ground survey station antenna under the ground fixation system calculated in the step (9)tECFTranslating the origin of coordinates from the geocentric to the antenna of the ground survey station, and calculating t1Position R of time satellite under the station center systemwCT：
RwCT＝MECF2CT*(RwECF-RtECF)
11. Direction-finding module for antenna of survey station
Calculating t from the position of the satellite under the satellite constellation calculated in the step (10)1The antenna pointing azimuth angle of the ground survey station at any moment: high and low angle
Figure BDA0002418106540000253
The method of the horizontal angle ψ is as follows:
high and low angles, horizontal angle. The high and low angles and the horizontal angle are defined under the station center system
Figure BDA0002418106540000254
To point to a vector RwCTAnd OCTXCTYCTAngle of plane, define RwCTVector sum of OCTZCTThe included angle is positive when the included angle is less than 90 degrees; the horizontal angle psi being the director vector RwCTAt OCTXCTYCTProjection of plane and OCTXCTAngle of axis, defined around OCTZCTShaft driven OCTXCTShaft clockwise steering pointing vector RwCTAt OCTXCTYCTThe projection of the surface is positive, as shown in fig. 5, and the antenna pointing azimuth is found according to this definition. Assuming that the position of the ground survey station is located at the origin of the station center system, the projection of the antenna pointing vector of the survey station under the station center system is RwCTRecord RwCTComprises the following steps:
Figure BDA0002418106540000255
calculation of elevation angle of antenna pointing azimuth
Figure BDA0002418106540000256
The horizontal angle ψ is:
Figure BDA0002418106540000261
Figure BDA0002418106540000262
and according to xCT、yCT、zCTPositive and negative of (2), angle of elevation
Figure BDA0002418106540000263
Dividing the horizontal angle psi into corresponding angle ranges to finish the forecasting process of the antenna pointing direction:
Figure BDA0002418106540000264
preferred embodiment 2:
the coordinate system required by the invention is as follows: the inertial system is a J2000.0 inertial coordinate system, and the earth fixation system is a WGS-84 coordinate system. The definition of the station center system is given below.
Standing heart system OCTXCTYCTZCT
The center of the station is defined as the origin OCTIs the ground antenna origin, the basic plane OCTXCTYCTThe surface is a local horizontal surface, OCTXCTPointing to true north, O, along the meridian of the local areaCTZCTVertical base plane pointing to zenith, OCTYCTDetermined by the right hand rule, as shown in fig. 4.
The calculation process of the present invention is detailed below:
the algorithm was simulated using MATLAB [ a ]0,e0,i0,Ω0,ω0,M0]Verification, the earth-related parameters and the station center are set as described above, and the ephemeris data of a certain type of satellite at the time 1, month 1, day 5 and day 4 of the UTC time 2019 are as follows:
Figure BDA0002418106540000265
Figure BDA0002418106540000271
six orbits of the satellite at 2019, 1, 6, 4 are recurred from 2019, 1, 5, 4, and the azimuth angle (elevation angle, horizontal angle) of the antenna pointing at 2019, 1, 6, 4 is calculated.
1. Input parameter singularity-free processing module
When the track eccentricity is very small (near-circular track, e ≈ 0), in order to avoid the occurrence of singular points in the calculation, the method is converted into 3 singular point-free variables, namely
ξ0＝e0 cos(ω0)
η0＝-e0 sin(ω0)
λ0＝ω0+M0
2. Input parameter normalization processing module
The normalization process is performed with respect to the recursion time dt,
dt＝t1-t0
Figure BDA0002418106540000272
for semi-major axis a0The normalization treatment is carried out, and the normalization treatment is carried out,
Figure BDA0002418106540000273
the normalized unit of the time is that,
Figure BDA0002418106540000274
wherein Ge is a gravitational constant;
the normalized unit of the length is Re-6378140 m
3. Perturbation item calculation module
Considering the perturbations of items J2 to J4 in the earth gravity field, including computing first order long term, first order short term, and second order long term terms, the expressions for the respective perturbation terms:
first order long term calculation:
Figure BDA0002418106540000281
Figure BDA0002418106540000282
Figure BDA0002418106540000283
the p is a radius, and the expression is p ═ a × (1-e)0 2) Mean angular velocity of track
Figure BDA0002418106540000284
J2＝1.624×10-3。
First order short period term calculation:
Figure BDA0002418106540000285
Figure BDA0002418106540000286
Figure BDA0002418106540000287
Figure BDA0002418106540000288
Figure BDA0002418106540000289
Figure BDA00024181065400002810
and u is calculated by a first-order long-term.
Figure BDA00024181065400002811
Figure BDA00024181065400002812
ξz1＝ξ0 cos(ωc1dtn)+η0 sin(ωc1dtn)
ηz1＝η0 cos(ωc1dtn)-ξ0 sin(ωc1dtn)
λz1＝λ0+(n+λc1)dtn
Figure BDA00024181065400002813
ωz1＝arc cos(ξz1/ez1)
Mz1＝mod(λz1-ωz1,2π)
Figure BDA0002418106540000291
u＝fz1+ωz1
Second-order long-term calculation:
Figure BDA0002418106540000292
Figure BDA0002418106540000293
Figure BDA0002418106540000294
said, J3＝2.5356×10-6,J4＝7.1022×10-6。
4. Recursion main formula calculation module
And for the introduced singularity-free variable, analyzing solution construction is carried out by combining a long-term item perturbation item and a short-period item perturbation item, and the long-period item perturbation is ignored. The track recursion main formula is as follows:
Figure BDA0002418106540000295
Figure BDA0002418106540000296
Figure BDA0002418106540000297
Figure BDA0002418106540000298
Figure BDA0002418106540000299
Figure BDA00024181065400002910
5. variable recovery module without singularity
After the calculation is finished, 3 singularity-free variables are restored
Figure BDA00024181065400002911
ωs＝arc cos(ξs/es)
Ms＝mod(λs-ωs,2π)
6. Normalized variable reduction module
Will be the long axis a of the satellite orbittReduction to conventional units:
as＝at×Re
a is describedsUnit: and m is selected.
The result obtained by the calculation in the steps is as follows:
Figure BDA0002418106540000301
7. inertial system satellite position module
Passing through t1Instantaneous number of satellite orbits [ a ] of times,es,is,Ωs,ωs,Ms]Calculating the component R of the position vector of the satellite in the J2000.0 coordinate systemwECIThe calculation method comprises the following steps:
RwECI＝Q*rp
wherein the rotation matrix Q is described in a 3-1-3 rotation order:
Figure BDA0002418106540000302
vector rp：
Figure BDA0002418106540000303
Wherein M is1True periapical angle:
Figure BDA0002418106540000304
the calculation result is as follows:
Figure BDA0002418106540000305
Figure BDA0002418106540000311
8. earth fixed satellite position module
According to a given target time t1And the position R of the satellite calculated in the step (7) in the inertial systemwECICalculating t1Position R of time satellite under earth's fixationwECFThe method comprises the following steps:
according to a given target time t1Calculating a second counting value t from epoch J2000.0 (1/12/2000) to a predetermined target timecInputting t1Year (year), month (month), day (day), hour (hour), minute (min), and second (sec) of time (UTC time), julian day JD is calculated:
Figure BDA0002418106540000312
wherein floor () is a round-down operation.
Calculating a second counting value t from epoch J2000.0 (1 month 1 day 12 in 2000) to a given target time according to the julian day JDc：
tc＝(JD-2455197.5)×86400+315547200
According to the second counting value t from the epoch J2000.0 (1 month, 1 day, 12 hours in 2000) to the given target timecAnd calculating a terrestrial rotation matrix ER, a nutation matrix NR and a time offset matrix PR, wherein the term is not considered in the invention because polar shift has little influence on the calculation of the conversion matrix. Calculating a transformation matrix M from an inertial system to a ground-based systemECI2ECF：
MECI2ECF＝ER*NR*PR
And counting according to the step (7)Calculated t1Position R of time satellite under inertial systemwECICalculating t1Position R of time satellite under earth's fixationwECF：
RwECF＝MECI2ECF*RwECI
The calculation result is as follows:
Figure BDA0002418106540000313
Figure BDA0002418106540000321
9. antenna position module of ground fixing system measuring station
The method is characterized in that the position R of the antenna of the ground measuring station under the ground fixation system is calculated according to the longitude, the latitude and the elevation of the given antenna of the ground measuring stationtECFThe method comprises the following steps:
and inputting longitude lon, latitude lat and elevation h of the antenna of the ground station.
Computing coordinate components G1, G2:
Figure BDA0002418106540000322
Figure BDA0002418106540000323
wherein f is the geometric oblateness of the earth ellipsoid, and f is 1/298.257.
Calculating the position R of the antenna of the ground survey station under the ground anchor systemtECF：
Figure BDA0002418106540000324
The calculation result is as follows:
Figure BDA0002418106540000325
10. satellite position module for station center system
T calculated according to the step (8)1Position R of time satellite under earth's fixationwECFAnd the position R of the ground survey station antenna calculated in the step (9) under the ground systemtECFCalculating t1Position R of time satellite under the system of the center of the stationwCTThe method comprises the following steps:
under the system of the station center, a conversion matrix M from the earth fixation system to the system of the station center is calculatedECF2CTDescribed as one rotation about the Z-axis of the earth fixation system and one rotation about the X-axis of the earth fixation system:
MECF2CT＝Rx(90°-lat)Rz(90°+lon)
wherein lon is the geographic longitude of the antenna of the ground survey station; lat is the geographic latitude of the ground station antenna.
Figure BDA0002418106540000331
Figure BDA0002418106540000332
And according to t calculated in the step (8)1Position R of time satellite under earth's fixationwECFThe position R of the ground survey station antenna under the ground fixation system calculated in the step (9)tECFTranslating the origin of coordinates from the geocenter to the antenna of the ground survey station, and calculating t1Position R of time satellite under the system of the center of the stationwCT：
RwCT＝MECF2CT*(RwECF-RtECF)
The calculation result is as follows:
Figure BDA0002418106540000333
11. direction-finding module for antenna of survey station
Calculating t from the position of the satellite under the satellite constellation calculated in the step (10)1The antenna pointing azimuth angle of the ground survey station at any moment: high and low angle
Figure BDA0002418106540000334
The method of the horizontal angle ψ is as follows:
high and low angles, horizontal angle. The high and low angles and the horizontal angle are defined under the station center system
Figure BDA0002418106540000335
To point to a vector RwCTAnd OCTXCTYCTAngle of plane, define RwCTVector sum of OCTZCTThe included angle is less than 90 degrees and is positive; the horizontal angle psi is the direction vector RwCTAt OCTXCTYCTProjection of plane and OCTXCTAngle of axis, defined around OCTZCTShaft driven by OCTXCTShaft clockwise steering pointing vector RwCTAt OCTXCTYCTThe projection of the surface is positive, as shown in fig. 5, and the antenna pointing azimuth is found according to this definition. Assuming that the position of the ground survey station is located at the origin of the station center system, the projection of the antenna pointing vector of the survey station under the station center system is RwCTRecord RwCTComprises the following steps:
Figure BDA0002418106540000341
calculation of elevation angle of antenna pointing azimuth
Figure BDA0002418106540000342
The horizontal angle ψ is:
Figure BDA0002418106540000343
Figure BDA0002418106540000344
and according to xCT、yCT、zCTPositive and negative of (2), angle of elevation
Figure BDA0002418106540000345
Dividing the horizontal angle psi into corresponding angle ranges to finish the forecasting process of the antenna pointing direction:
Figure BDA0002418106540000346
the calculation result is as follows:
Figure BDA0002418106540000347
in the description of the present application, it is to be understood that the terms "upper", "lower", "front", "rear", "left", "right", "vertical", "horizontal", "top", "bottom", "inner", "outer", and the like indicate orientations or positional relationships based on those shown in the drawings, merely for convenience of description and simplicity of description, and do not indicate or imply that the device or element being referred to must have a particular orientation, be constructed in a particular orientation, and be operated, and therefore, are not to be construed as limiting the present application.
Those skilled in the art will appreciate that, in addition to implementing the systems, apparatus, and various modules thereof provided by the present invention in purely computer readable program code, the same procedures can be implemented entirely by logically programming method steps such that the systems, apparatus, and various modules thereof are provided in the form of logic gates, switches, application specific integrated circuits, programmable logic controllers, embedded microcontrollers and the like. Therefore, the system, the device and the modules thereof provided by the present invention can be considered as a hardware component, and the modules included in the system, the device and the modules thereof for implementing various programs can also be considered as structures in the hardware component; modules for performing various functions may also be considered to be both software programs for performing the methods and structures within hardware components.
The foregoing description of specific embodiments of the present invention has been presented. It is to be understood that the present invention is not limited to the specific embodiments described above, and that various changes or modifications may be made by one skilled in the art within the scope of the appended claims without departing from the spirit of the invention. The embodiments and features of the embodiments of the present application may be combined with each other arbitrarily without conflict.

Claims (9)

1. A method for calculating an azimuth angle of a ground measurement and control antenna pointing to a satellite is characterized by comprising the following steps:
using ephemeris time t0、t0Number of orbital flat at time [ a ]0,e0,i0,Ω0,ω0,M0]And the geographic position information of the ground survey station antenna is finally converted into t through the correlation calculation of the satellite orbit and the conversion calculation of a plurality of correlation coordinate systems1Pointing azimuth angle of antenna under standing center system at time: high and low angle
Figure FDA0003404198430000011
The horizontal angle psi is used for completing the prediction process of the satellite pointing direction by the ground station-measuring antenna;
the method comprises the following steps:
the input parameter non-singular processing steps:
when the track eccentricity is small, i.e. near circular track, e ≈ 0, to avoid the occurrence of singularities in the calculation, we turn to 3 singularity-free variables, i.e.:
ξ0＝e0 cos(ω0)
η0＝-e0 sin(ω0)
λ0＝ω0+M0
input parameter normalization processing:
normalization processing is performed for the recursion time dt:
dt＝t1-t0
Figure FDA0003404198430000012
wherein the content of the first and second substances,
t0recursion of the starting moment for the track;
t1representing the track recursion end time;
dtnthe result of normalization processing of the recursion time is shown, and the subscript n shows normalization;
for semi-major axis a0And (3) carrying out normalization treatment:
Figure FDA0003404198430000013
the normalized unit of the time is
Figure FDA0003404198430000021
Wherein the content of the first and second substances,
ge is a gravitational constant;
re represents the Earth's equatorial radius;
perturbation item calculation step:
considering the perturbations of items J2-J4 in the earth gravity field, including the calculation of first-order long-term, first-order short-term, and second-order long-term terms, the expressions for the individual perturbation terms are as follows:
first order long term calculation:
Figure FDA0003404198430000022
Figure FDA0003404198430000023
Figure FDA0003404198430000024
wherein the content of the first and second substances,
Ω1representing the first-order long period variable quantity of the right ascension of the satellite orbit intersection point;
ω1representing the first-order long period variation of the argument of the satellite orbit in the near place;
λ1representing the first-order long period variation of the introduced singularity-free variable λ is calculated;
p is a radius, and the expression is p ═ a × (1-e)0 2) Mean angular velocity of track
Figure FDA0003404198430000025
J2＝1.624×10-3；
First order short period term calculation:
Figure FDA0003404198430000026
Figure FDA0003404198430000027
Figure FDA0003404198430000028
Figure FDA0003404198430000029
Figure FDA00034041984300000210
Figure FDA00034041984300000211
wherein the content of the first and second substances,
Figure FDA0003404198430000031
representing the first-order short period variable quantity of the semi-major axis of the satellite orbit;
Figure FDA0003404198430000032
representing the first-order short-period variation of the satellite orbit inclination angle;
Figure FDA0003404198430000033
the first-order short-period variation of the right ascension of the satellite orbit intersection point is represented;
Figure FDA0003404198430000034
respectively representing the first-order short-period variable quantity of three singularity-free variables xi, eta and lambda introduced by calculation;
the u is calculated by a first-order long-term, and the calculation process is as follows:
ξz1＝ξ0cos(ω1dtn)+η0sin(ω1dtn)
ηz1＝η0cos(ω1dtn)-ξ0sin(ω1dtn)
λz1＝λ0+(n+λ1)dtn
Figure FDA0003404198430000035
ωz1＝arc cos(ξz1/ez1)
Mz1＝mod(λz1-ωz1，2π)
Figure FDA0003404198430000036
u＝fz1+ωz1
wherein the content of the first and second substances,
ξz1、ηz1、λz1respectively representing the integral results of three singularity-free variables xi, eta and lambda introduced by calculation according to the first-order long-term variation;
ez1an integral result which represents the satellite orbit eccentricity and takes the first-order long period variable quantity as an integral quantity;
ωz1an integral result which represents the argument of the satellite orbit near place and takes the first-order long period variable quantity as integral quantity;
Mz1an integration result which represents the mean-near point angle of the satellite orbit and takes the first-order long period variable quantity as an integral quantity;
fz1an integral result which represents the true near point angle of the satellite orbit and takes the first-order long period variable quantity as an integral quantity;
second-order long-term calculation:
Figure FDA0003404198430000041
Figure FDA0003404198430000042
Figure FDA0003404198430000043
wherein the content of the first and second substances,
Ω2representing the second-order long-period variation of the right ascension of the satellite orbit;
ξ2、λ2respectively representing the second-order long-period variable quantity of the introduced singularity-free variables xi and lambda;
J3＝2.5356×10-6,J4＝7.1022×10-6。
2. the method of claim 1, further comprising:
and (3) calculating a recursion main formula:
for introduced singularity-free variables, an analytic solution structure is carried out by combining a long-term item perturbation item and a short-period item perturbation item, the long-period item perturbation is ignored, and a track recursion main formula is as follows:
Figure FDA0003404198430000044
Figure FDA0003404198430000045
Figure FDA0003404198430000046
Figure FDA0003404198430000047
Figure FDA0003404198430000048
Figure FDA0003404198430000049
wherein the content of the first and second substances,
αtrepresents t1At the moment, normalizing the result of the satellite orbit semi-major axis;
isrepresents t1At that time, the satellite orbit inclination;
Ωsdenotes t1At the moment, the rising point of the satellite orbit is right ascension;
ξs、ηs、λsrespectively represent t1At the moment, the values of the three singularity-free variables xi, η, λ introduced are calculated.
3. The method of claim 2, further comprising:
a singular point-free variable reduction step:
after the calculation is finished, 3 singularity-free variables are restored
Figure FDA0003404198430000051
ωs＝arc cos(ξs/es)
Ms＝mod(λs-ωs,2π)。
4. The method of claim 3, further comprising:
and (3) normalization variable reduction step:
will be the long axis a of the satellite orbittReduction to conventional units:
as＝at×Re
a is describedsUnit: and m is selected.
5. The method of claim 4, further comprising:
calculating the position of the inertial system satellite:
input t1Instantaneous number of satellite orbits [ a ] of times,es,is,Ωs,ωs,Ms]Outputting the component R of the position vector of the satellite in the J2000.0 coordinate systemwECIThe calculation method comprises the following steps:
RwECI＝Q*rp
wherein the rotation matrix Q is described in a 3-1-3 rotation order:
Figure FDA0003404198430000052
vector rp：
Figure FDA0003404198430000053
Wherein M is1True proximal angle:
Figure FDA0003404198430000054
。
6. the method of claim 5, further comprising:
and a step of calculating the position of the earth-fixed satellite:
inputting a target time t1And said calculated position R of the satellite in the inertial systemwECIOutput t1Position R of time satellite under earth's fixationwECFThe specific method comprises the following steps:
according to a given target time t1Calculating the second counting value t from epoch J2000.0 (1/12/2000) to the given target timecInputting t1Year, month, day, hour, minute, second of the moment, calculate julian day JD:
Figure FDA0003404198430000061
wherein, floor () is a round-down operation;
calculating a second count value t from epoch J2000.0 to a given target time based on the julian day JDc：
tc＝(JD-2455197.5)×86400+315547200
According to the calculated epoch J2000.0 to the second counting value t of the given target timecCalculating a rotation matrix ER, a nutation matrix NR and a precision matrix PR of the earth, and calculating a conversion matrix M from an inertia system to a ground-fixed systemECI2ECF：
MECI2ECF＝ER*NR*PR
And according to t obtained by the calculation1Position R of time satellite under inertial systemwECICalculating t1Position R of time satellite under earth's fixationwECF：
RwECF＝MECI2ECF*RwECI。
7. The method of claim 6, further comprising:
and the antenna position calculation step of the earth-fixed system survey station:
according to the longitude, the latitude and the elevation of the given ground survey station antenna, the position R of the ground survey station antenna under the ground fixation system is calculatedtECFThe specific method comprises the following steps:
inputting longitude lon, latitude lat and elevation h of an antenna of the ground station;
computing coordinate components G1, G2:
Figure FDA0003404198430000071
Figure FDA0003404198430000072
wherein f is the geometric oblateness of the earth ellipsoid, and f is 1/298.257;
calculating the position R of the ground survey station antenna under the ground fixing systemtECF：
Figure FDA0003404198430000073
。
8. The method of claim 7, further comprising:
and a step of calculating the position of the satellite in the station center system:
inputting calculated t1Position R of time satellite under earth's fixationwECFAnd said calculated position R of the ground station antenna under the ground anchor systemtECFOutputting the position R of the satellite under the station center systemwCTThe specific method comprises the following steps:
under the station-centric system, a conversion matrix M from the earth-fixed system to the station-centric system is calculatedECF2CTDescribed as one rotation about the Z axis of the earth fixation and one rotation about the X axis of the earth fixation:
MECF2CT＝Rx(90°-lat)Rz(90°+lon)
wherein lon is the geographic longitude of the antenna of the ground survey station; lat is the geographical latitude of the ground survey station antenna;
Figure FDA0003404198430000074
Figure FDA0003404198430000075
and according to said t1Position R of time satellite under earth fixed systemwECFPosition R of the ground station antenna under the ground anchortECFTranslating the origin of coordinates from the geocentric to the antenna of the ground survey station, and calculating t1Position R of time satellite under the system of the center of the stationwCT：
RwCT＝MECF2CT*(RwECF-RtECF)。
9. The method of claim 8, further comprising:
the method comprises the following steps of:
inputting the position of the satellite under the station center system and outputting t1The antenna pointing azimuth angle of the ground survey station at any moment: high and low angle
Figure FDA00034041984300000810
The horizontal angle psi is specifically as follows:
the high and low angles and the horizontal angle are defined under the station center system OCTXCTYCTZCTIn, OCTRepresenting the origin of the coordinate system, XCTX-axis, Y, representing a coordinate systemCTY-axis, O, representing a coordinate systemCTXCTYCTThe plane, elevation angle, of X-axis and Y-axis of the coordinate system
Figure FDA0003404198430000081
To point to a vector RwCTAnd OCTXCTYCTAngle of plane, define RwCTVector sum of OCTZCTThe included angle is less than 90 degrees and is positive; the horizontal angle psi being the director vector RwCTAt OCTXCTYCTProjection of plane and OCTXCTAngle of axis, defined around OCTZCTShaft driven by OCTXCTShaft clockwise steering pointing vector RwCTAt OCTXCTYCTThe projection of the surface is positive, the antenna direction azimuth angle is obtained according to the definition, and the projection of the antenna direction vector of the measuring station under the station center system is R if the position of the ground measuring station is positioned at the origin of the station center systemwCTRecord RwCTComprises the following steps:
Figure FDA0003404198430000082
wherein the content of the first and second substances,
xCT、yCT、zCTrespectively representing the pointing vectors R of antennas of the stationswCTAt the center of the standOCTXCTYCTZCTCoordinate components of corresponding X-axis, Y-axis and Z-axis;
calculation of elevation angle of antenna pointing azimuth
Figure FDA0003404198430000083
The horizontal angle ψ is:
Figure FDA0003404198430000084
Figure FDA0003404198430000085
and according to xCT、yCT、zCTPositive and negative of (2), angle of elevation
Figure FDA0003404198430000086
Dividing the horizontal angle psi into corresponding angle ranges to finish the process of forecasting the antenna pointing direction:
if z isCTMore than or equal to 0, the height and angle will be changed
Figure FDA0003404198430000087
Division into ranges
Figure FDA0003404198430000088
If xCT≥0,yCTThe horizontal angle psi is divided into ranges
Figure FDA0003404198430000089
If xCT＜0,yCTThe horizontal angle psi is divided into ranges
Figure FDA0003404198430000091
If xCT＜0,yCT< 0, divide the horizontal angle psi into the range
Figure FDA0003404198430000092
If xCT≥0,yCTIf < 0, the horizontal angle psi is divided into ranges
Figure FDA0003404198430000093
CN202010197399.2A 2020-03-19 2020-03-19 Azimuth angle calculation method for ground measurement and control antenna pointing satellite Active CN111427002B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010197399.2A CN111427002B (en) 2020-03-19 2020-03-19 Azimuth angle calculation method for ground measurement and control antenna pointing satellite

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010197399.2A CN111427002B (en) 2020-03-19 2020-03-19 Azimuth angle calculation method for ground measurement and control antenna pointing satellite

Publications (2)

Publication Number Publication Date
CN111427002A CN111427002A (en) 2020-07-17
CN111427002B true CN111427002B (en) 2022-07-12

Family

ID=71548179

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010197399.2A Active CN111427002B (en) 2020-03-19 2020-03-19 Azimuth angle calculation method for ground measurement and control antenna pointing satellite

Country Status (1)

Country Link
CN (1) CN111427002B (en)

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112329202B (en) * 2020-09-30 2023-08-15 北京空间飞行器总体设计部 Optimization implementation method of antenna pointing algorithm of circulator by Mars
CN112948741B (en) * 2021-02-04 2023-02-28 上海卫星工程研究所 Method and system for calculating visible arc section of deep space probe
CN114137573A (en) * 2021-09-16 2022-03-04 北京微纳星空科技有限公司 Satellite measurement and control station, satellite measurement and control method, equipment and storage medium
CN114002713A (en) * 2021-10-20 2022-02-01 上海航天空间技术有限公司 Satellite orbit parameter recursion processing and forecasting system
CN114002710A (en) * 2021-10-20 2022-02-01 上海航天空间技术有限公司 On-satellite orbit position autonomous prediction method for small-eccentricity low-orbit satellite
CN114826438A (en) * 2022-03-30 2022-07-29 中国空间技术研究院 Method for determining installation position of ground measurement and control antenna
CN116659543A (en) * 2023-06-21 2023-08-29 中国人民解放军61540部队 Satellite position and attitude estimation method and device based on remote sensing satellite orbit number

Family Cites Families (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPS61502740A (en) * 1984-07-20 1986-11-27 ヒユ−ズ・エアクラフト・カンパニ− Rotating stabilized satellite with nutation control subsystem
CN102230969B (en) * 2011-03-22 2013-05-01 航天恒星科技有限公司 Long-time independent maintenance method of inter-satellite link in satellite constellation
CN104729457B (en) * 2015-04-16 2017-04-12 哈尔滨工业大学 Method for determining position of sun relative to near-earth microsatellites
CN105158777B (en) * 2015-07-31 2017-08-29 上海卫星工程研究所 The data source generation method positioned for passive direction finding
CN107656293A (en) * 2017-08-30 2018-02-02 清华大学 The 14 parameter broadcast ephemeris satellite position methods of estimation based on Nonsingular orbital elements
CN108974395B (en) * 2018-06-21 2019-11-15 中国人民解放军战略支援部队航天工程大学 Extraterrestrial target based on the driving of sky-based laser platform becomes rail calculation method and its device
CN110838864A (en) * 2018-08-19 2020-02-25 南京理工大学 Unattended satellite ground station tracking control system
CN110816896B (en) * 2019-10-28 2021-06-11 中国空间技术研究院 Satellite on-satellite simple orbit extrapolation method

Also Published As

Publication number Publication date
CN111427002A (en) 2020-07-17

Similar Documents

Publication Publication Date Title
CN111427002B (en) Azimuth angle calculation method for ground measurement and control antenna pointing satellite
US5957982A (en) Method and system for space navigation
CN112629543A (en) Orbit planning method for large elliptical orbit and small-inclination-angle circular orbit
CN111427003A (en) Pointing guidance system of ground survey station antenna to satellite
CN112713922A (en) Visibility rapid forecasting algorithm of multi-beam communication satellite
CN112649006A (en) Orbit planning method for sun synchronous circular orbit
CN111427004A (en) Coordinate conversion method suitable for pointing of ground survey station antenna to satellite
CN112014869A (en) Astronomical navigation-based inter-satellite link autonomous navigation method and system
CN111427001A (en) Target positioning method suitable for pointing satellite by ground survey station antenna
CN111679242A (en) Ground antenna guiding method suitable for pointing to in-orbit spacecraft
Liu et al. Autonomous orbit determination and timekeeping in lunar distant retrograde orbits by observing X‐ray pulsars
Barbee et al. Guidance and navigation for rendezvous and proximity operations with a non-cooperative spacecraft at geosynchronous orbit
Liu et al. Guidance and control technology of spacecraft on elliptical orbit
Gustavsson Development of a MatLab based GPS constellation simulation for navigation algorithm developments
CN112013834B (en) Astronomical navigation-based inter-satellite link autonomous recovery method and system
CN111366126A (en) System for calculating apparent vector of satellite pointing by ground survey station antenna
Lincoln et al. Components of a vision assisted constrained autonomous satellite formation flying control system
CN111427000A (en) Target view vector determination method suitable for pointing of ground survey station antenna to satellite
CN114002710A (en) On-satellite orbit position autonomous prediction method for small-eccentricity low-orbit satellite
Wang et al. Autonomous orbit determination using pulsars and inter-satellite ranging for Mars orbiters
CN112394381A (en) Full-autonomous lunar navigation and data communication method based on spherical satellite
Parker et al. Navigation between geosynchronous and lunar L1 orbiters
Tong et al. Relative motion control for autonomous rendezvous based on classical orbital element differences
Zhou A study for orbit representation and simplified orbit determination methods
Peng et al. UPF based autonomous navigation scheme for deep space probe

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant