CN106570270B - A kind of more star combined covering characteristic fast determination methods of System of System Oriented design - Google Patents
A kind of more star combined covering characteristic fast determination methods of System of System Oriented design Download PDFInfo
- Publication number
- CN106570270B CN106570270B CN201610972224.8A CN201610972224A CN106570270B CN 106570270 B CN106570270 B CN 106570270B CN 201610972224 A CN201610972224 A CN 201610972224A CN 106570270 B CN106570270 B CN 106570270B
- Authority
- CN
- China
- Prior art keywords
- satellite
- equator
- coverage
- point
- intersection
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 56
- 238000013461 design Methods 0.000 title claims abstract description 20
- 238000005070 sampling Methods 0.000 claims description 102
- 238000004364 calculation method Methods 0.000 claims description 49
- 238000012935 Averaging Methods 0.000 claims description 3
- 238000004088 simulation Methods 0.000 description 13
- 238000011156 evaluation Methods 0.000 description 5
- 238000005457 optimization Methods 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- 238000006243 chemical reaction Methods 0.000 description 2
- 238000010276 construction Methods 0.000 description 2
- 230000009897 systematic effect Effects 0.000 description 2
- 238000013524 data verification Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000003203 everyday effect Effects 0.000 description 1
- 230000000737 periodic effect Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Navigation (AREA)
Abstract
The present invention relates to a kind of more star combined covering characteristic fast determination methods of System of System Oriented design, and the nodal period of satellite and the equatorial plane is determined using three track characteristic parameters；Difference of longitude on determining satellite adjacent orbit sub-satellite track under the line；Determine all intersection point longitudes of Track of Sub-Satellite Point and equator in one day；Equator sampled point longitude is obtained, determines that satellite load is to the coverage condition of equator sampled point at equatorial node；And using all intersection point longitudes for the initial position and satellite same day sub-satellite track and equator for determining satellite next day, then satellite load is obtained in one day to the coverage condition of equator sampled point；Obtain each star equatorial node longitude and to equator sampled point coverage condition；Judge whether equator sampled point covers completely, obtains more star combined covering time index values.The present invention can greatly reduce orbit prediction and observation coverage amount, and have certain accuracy guarantee, while calculate simply, and Project Realization is easy.
Description
Technical Field
The invention relates to a method for quickly calculating the global coverage time of a multi-satellite combination, which is particularly suitable for quickly finishing the simulation calculation of the coverage time of a plurality of satellite combination modes under the conditions of planning design, efficiency evaluation and the like of a complex spacecraft system so as to support the iterative optimization design of a system scheme.
Background
The multi-satellite combination coverage characteristics comprise various characteristic indexes such as multi-satellite combination global or designated area coverage time, multi-satellite combination designated area or designated place revisit time and the like. Generally, for system design simulation, the global coverage time of multi-satellite combination is one of the most representative characteristic indexes.
With the vigorous development of the aerospace field, the aerospace remote sensing application requirements and the number and types of remote sensing satellites are rapidly increased, one-satellite multi-use and multi-satellite public use become the basis standards for systematic design requirement input and application efficiency evaluation, and the systematic design construction of space-based resources and the coordinated sharing of resources become the important trends of aerospace system planning, construction and application, so that the problems of spacecraft system optimization design and application evaluation are generated. The biggest difficulty of the problem is that for system planning design, the required number is huge, the number and load types of the satellites forming the system are far higher than those of the situation that a plurality of satellites form a constellation under the general condition, multiple satellite load combination application modes exist according to specific application requirements, when the number of required items is increased and observation tasks conflict, the satellites need to be combined in multiple ways to meet the requirements as much as possible, and the traversal calculation amount is large. Meanwhile, considering that six orbits are empirical rough design in the system design stage, especially the observation capability of the on-orbit satellite needs to be comprehensively considered, long-time orbit prediction can be carried out, and the uncertainty of remote sensing task design in application can influence the coverage numerical simulation precision. Therefore, accurate numerical simulation calculation does not meet the calculation requirement of the situation, but a simplified calculation method which has good universality and can obviously improve the coverage simulation efficiency is needed.
Currently, a numerical simulation method is adopted for solving the problem of single-star or multi-star coverage characteristic calculation, and is generally used for simulating and evaluating multi-star coverage performance simulation. The method comprises the following specific steps: (1) rasterizing a coverage area based on the coverage calculation precision requirement; (2) setting an orbit forecasting simulation step length, carrying out orbit forecasting on each satellite according to all satellite orbit characteristic parameters (six orbits) participating in coverage calculation, and calculating a sub-satellite point track; (3) and calculating coverage characteristic indexes such as revisit, coverage time and the like by combining information such as the load width, observation angle and the like of each satellite. The calculation method has the advantages that the calculation result is more accurate under the conditions that the grid granularity is fine, the parameters such as the track, the breadth, the observation angle and the like are accurate, and the perturbation model is reasonably selected; the disadvantages are also obvious, namely, the calculation amount is large, the time consumption is long, and especially in the case of a large number of satellites participating in calculation or a large number of satellite combination types aiming at application, the simulation calculation amount of the system design is too large.
In addition, in engineering, a common coverage characteristic experience estimation method is provided for a plurality of satellites of the same type in the same orbital plane, and the method is commonly used for estimating the number of satellites of the same type in the same orbit meeting certain coverage requirements in a territorial scope and designing phase distribution. The method comprises the following specific steps: (1) calculating the number of turns around the earth or the orbit period every day according to the orbit parameters, and estimating the effective observation orbit number of the coverage area according to the orbit type (such as morning and afternoon stars); (2) then, the required total track number and the phase difference of adjacent tracks are calculated approximately according to the load width, the boundary range of the observation area and the requirement on the revisit time of the area; (3) and calculating the phase distribution and the coverage characteristic parameters of the multiple stars by combining the load width of each star. The approximate calculation method assumes ideal conditions, is a calculation situation of roughly calculating a specific region in a short time, cannot be quickly converted into different tracks and different loads, and has insufficient universality and insurable precision.
In system efficiency evaluation, multi-star coverage characteristic comparison analysis is a basic index calculation problem. When the calculated amount is not large and no requirement is required for calculation time consumption, the first numerical calculation method mentioned above is adopted, so that the requirement can be met; however, when the number of satellites involved in calculation is large, the time span spanned by system design is long, and the required number to be evaluated is large, except for assuming that an ideal satellite combination with the same orbital plane can adopt a second quick calculation method, for the conditions of different orbits and different load widths, if the system performance indexes with the coverage characteristics of one-satellite multi-use, multi-satellite public use and the like need to be evaluated and calculated, the simulation time is lengthened due to the numerical calculation complexity; meanwhile, in the planning and design stage, the orbit parameter design of the satellite generally comprises the steps of firstly setting an empirical value and then carrying out optimization and adjustment according to system simulation, so that the calculation still needs repeated iteration, and the corresponding calculation needs longer time.
For the whole system, when the areas to be calculated are distributed dispersedly in the global range and the number is variable, it is difficult to quickly estimate the coverage of the areas by different tracks and different load widths, so that it is necessary to improve the method, make conditional assumptions for different situations, and design a corresponding simplified calculation method and a general judgment criterion.
Disclosure of Invention
The technical problem to be solved by the invention is as follows: the method overcomes the defects of the prior art, can greatly reduce the calculation amount of track forecasting and observation coverage, has certain precision guarantee, and is simple in calculation and easy in engineering realization.
The technical solution of the invention is as follows: a system design-oriented multi-satellite combination coverage characteristic rapid determination method is characterized by comprising the following steps:
(1) determining the intersection point period of the satellite and the equatorial plane by utilizing three orbit characteristic parameters of an orbit semi-major axis, an orbit eccentricity and an orbit inclination angle;
(2) determining the longitude difference of the track of the subsatellite point of the adjacent orbit of the satellite on the equator by using the earth rotation angular velocity and the intersection point cycle obtained in the step (1);
(3) taking the intersection points of the satellite initial orbit subsatellite points and the equator as the initial positions of the satellites, and determining all the longitude of the intersection points of the satellite subsatellite point tracks and the equator in one day by the intersection point period of the satellites and the equator plane and the longitude difference of the satellite subsatellite point tracks of the adjacent orbits in the step (2) on the equator;
(4) setting sampling granularity under the full coverage condition of the equator, obtaining longitude of sampling points of the equator, and determining the coverage condition of satellite loads at the intersection points of the equator on the sampling points of the equator by using the satellite load breadth and all obtained longitude of intersection points of the equator;
(5) taking the integral multiple of the satellite intersection period and taking the number of the satellite intersection period as a unit of day, determining the initial position of the next day of the satellite and all the intersections longitude of the orbit of the satellite under the current day and the equator by using the steps (1) to (3), and obtaining the coverage condition of the satellite load on the equator sampling point in one day by using the satellite load width and the equator sampling point longitude in the step (4);
(6) repeating the steps (1) to (5) by utilizing three orbit characteristic parameters of a multi-star semi-major axis, orbit eccentricity and orbit inclination to obtain the longitude of each star equator intersection point and the coverage condition of each equator sampling point;
(7) and (5) judging whether the equator sampling points are completely covered or not, if not, repeating the step (5), and if the equator sampling points are completely covered and the days are completely recorded, obtaining the multi-satellite combination coverage time index value.
The method comprises the following steps that (1) the intersection point period of the satellite and the equatorial plane is determined by utilizing three orbit characteristic parameters including orbit semimajor axis, orbit eccentricity and orbit inclination, and the intersection point period is obtained by the following method:
the period of the intersection points is set to
Wherein:
where n is the satellite motion averaging period, i.e.
In the formula,is the argument of the near placeThe position is radian;is a flat proximal angle, with units of radians; a is a semi-major axis of the track, and the unit is meter; e is track eccentricity, i is track inclination, J21.082625829965505e-3, second-order band harmonic coefficients of earth's gravitational potential, Re6378137.0, the radius of the earth's equator, in meters; mu.sE3.98600436e +14, is the gravitational constant of the earth, in meters3Second/second2。
The method for determining the longitude difference of the track of the subsatellite point of the adjacent orbit of the satellite in the step (2) on the equator comprises the following steps:
let satellite xjDifference in longitude Δ λ between adjacent tracks on the equatorj，
Wherein,
the rate of change of the right ascension at the intersection point;is the rotational angular velocity of the earth; a is a semi-major axis of the track, and the unit is meter; e is track eccentricity and i is track inclination.
In the step (3)
(1) The method for determining the initial position of each satellite comprises the following steps:
let the initial position of the jth star be xloc(j)，xloc(j) Is a random number between (0,2 pi) and is greater than or equal to the minimum load width W in the multi-satellite combinationminExpressed in radians; obtaining the initial position of each satellite in the same way;
(2) the longitude determination method for all intersection points of the satellite subsatellite point track and the equator in one day comprises the following steps:
Xj＝{xloc(j)+(n-1)·Δλjn belongs to N and N is less than or equal to cjIn which xloc(j) Is the initial position of the jth satellite, n is the ordinal number of the intersection of the orbit of the satellite's subsatellite point and the equator, cjIs a satellite xjThe number of the intersection periods in one day,meaning that the rounding is done down,the period of the intersection of the satellite obtained in step (1) with the equatorial plane is expressed in hours.
In the step (4)
(1) The longitude coverage condition of the equatorial sampling points is expressed by the following method:
carrying out rasterization equipartition sampling on the equator, wherein under the condition of N equipartition, N +1 sampling points exist, and the position of the p-th sampling point is recorded as NpTaking the coverage of the satellite width on all sampling points as a coverage condition judgment basis, and recording the number of uncovered sampling points as NIs prepared fromAnd the threshold value of the number of uncovered sampling points for finishing the calculation of the full coverage cycle is recorded as NtolIntroducing a flag function flag (p) at each sampling point,
(2) the method for determining the coverage condition of the satellite load at the equator intersection point to the equator sampling point comprises the following steps:
with satellite xjFor the coverage of the c-th intersection on day d as an example, first, the satellite x is calculatedjCurrent coverage interval ofWhereinxloc(j) Is the initial position of the jth star in radians; delta lambdajIs a satellite xjLongitude differences of adjacent tracks on the equator are in radian units, and whether N +1 equator sampling points are covered or not is sequentially judged;
the step of whether the equator sampling point is covered or not is as follows: if the current sampling point NpIf the current sampling point is not covered, namely the flag function value flag (p) of the current sampling point is 0, the sampling point N is judgedpWhether or not in a coverage areaWithin the range: if N is presentpIn the interval, the flag function value flag (p) of the current sampling point is changed to 1, and the number of uncovered sampling points is recorded as NIs prepared from＝NIs prepared from-1; if N is presentpAnd (4) not in the interval, performing no operation, circulating for N +1 times, and finishing the coverage condition determination of each equatorial sampling point.
The method for determining the coverage condition of the satellite load on the equator sampling point in one day in the step (5) comprises the following steps:
for satellite xjAnd (4) determining the coverage of c intersection points on the day d by using the method for determining the coverage of the satellite load on the equator sampling point at the equator intersection point in the step (4) And (4) performing coverage judgment on the equator sampling point in the step (4) in each coverage interval in a circulating manner to obtain the coverage condition of the satellite load on the equator sampling point in one day.
The method for determining the longitude and the coverage of each star equatorial sampling point in the step (6) comprises the following steps:
for the coverage condition of the day d, firstly, the initial position (unit is radian) of each star and the longitude difference of the adjacent track of each star on the equator, which are required by the calculation of the step (4) and the step (5), are obtained by using the steps (1) to (3), and then, the step (5) is repeatedly used for each star, namely, the longitude of each star equator sampling point and the coverage condition of the equator sampling point on the day d are obtained.
The method for determining the number of days d covered by the multi-satellite combination coverage time index value in the step (7) comprises the following steps:
assuming that the initial state d is 0, after step (6) is completed, it is determined whether to cover all the equatorial sampling points, i.e. NIs prepared fromWhether or not less than NtolIn which N isIs prepared fromNumber of uncovered samples, NtolThreshold number of uncovered samples to end the full coverage cycle, if NIs prepared fromGreater than NtolIf d is d +1, the step (6) is executed in a loop, and if N is not equal to d +1Is prepared fromIs less than or equal to NtolAnd obtaining the number of days d covered by the multi-satellite combined covering time index value.
Compared with the prior art, the invention has the advantages that:
(1) the method simplifies the global satellite down-pointing forecast calculation of the satellite orbit by utilizing the periodic characteristics of the earth observation satellite operation and the calculation method of the intersection point longitude of the satellite operation orbit and the equator, thereby greatly reducing the orbit forecast calculation amount of the service global coverage calculation; by using the intersection point period calculation model considering J2 perturbation, the orbit intersection point longitude is ensured to have higher calculation precision, and the application requirement of system design simulation can be met.
(2) According to the method, the full coverage characteristic calculation of the satellite load on the equator is used for replacing the judgment calculation of the load global coverage geographic grid, namely, the one-dimensional coverage condition calculation of the width on the equator circumference is used for replacing the two-dimensional calculation of the global coverage, so that the calculation dimensionality of the circular judgment of the geographic grid is reduced, and the multi-step coordinate system conversion calculation of the two-dimensional coordinate from the equator to the two-polar region sampling point in the grid calculation is avoided; compared with the common two-dimensional sampling point coverage judgment, the method avoids the calculation error caused by the coordinate system conversion, has simpler and more convenient calculation and is easy to realize engineering.
(3) The intersection point period information obtained by the invention can be used as the basis for further calculating the coverage characteristic calculation of the area. When the coverage characteristic analysis of the system-oriented satellite load specific area is carried out, the lowest dimensional circle of the earth where the area is located can be used as the equator circle in the determination method, the longitude difference between adjacent orbits on the corresponding dimensional circle is obtained by calculating the period of intersection points and combining the rotation angular velocity and the orbit inclination angle of the earth, and the dimensional circle is sampled by utilizing the subsequent steps, so that the full coverage characteristic of the area can be rapidly determined, the judgment and calculation of the area coverage grids can be reduced, and the engineering application is facilitated.
In a word, the method passes the actual data verification in the spacecraft system simulation and performance evaluation system, is feasible, and is easy to realize by engineering technology, so that the method has practicability.
Drawings
Fig. 1 is a flow chart of fast determining the coverage characteristics of the multi-satellite combination based on the system-oriented design.
Detailed Description
Example 1:
as shown in fig. 1, the method comprises the following specific steps:
(1) and determining the intersection point period of the satellite and the equatorial plane by using three orbit characteristic parameters of the orbit semimajor axis, the orbit eccentricity and the orbit inclination angle.
The period of the intersection points is set to
Wherein:
where n is the satellite motion averaging period, i.e.
In the formula,is the argument of the near location (in radians),is mean paraxial angle (in radians), a is the semi-major axis of the track (in meters), e is the eccentricity of the track, i is the inclination of the track, J21.082625829965505e-3, second-order band harmonic coefficients of earth's gravitational potential, Re6378137.0, the equatorial radius of the earth (in meters), μE3.98600436e +14, is the gravitational constant of the earth (unit is meter)3Second/second2)。
(2) And (3) determining the longitude difference of the track of the subsatellite point of the adjacent orbits of the satellite on the equator by using the rotational angular velocity of the earth and the intersection point period obtained in the step (1).
Let satellite xjDifference in longitude of adjacent tracks on equatorΔλj，
Wherein,
the rate of change of the right ascension at the intersection point;is the rotational angular velocity of the earth; a is the track semimajor axis (in meters), e is the track eccentricity, and i is the track inclination.
(3) And (3) taking the intersection points of the initial orbital sub-satellite points of the satellite and the equator as the initial positions of the satellites, and determining all the longitude of the intersection points of the orbital sub-satellite tracks of the satellite and the equator in one day by the intersection period of the satellites and the equatorial plane and the longitude difference of the orbits of the adjacent orbits of the satellite in the step (2) on the equator.
Firstly, the method for determining the initial position of each satellite comprises the following steps:
let the initial position of the jth star be xloc(j)，xloc(j) Is a random number between (0,2 pi) and is greater than or equal to the minimum load width W in the multi-satellite combinationminExpressed in radians; and obtaining the initial position of each star in the same way.
Then, the longitude determination method of all intersection points of the satellite subsatellite point track and the equator in one day comprises the following steps:
Xj＝{xloc(j)+(n-1)·Δλjn belongs to N and N is less than or equal to cjIn which xloc(j) Is the initial position of the jth satellite, n is the ordinal number of the intersection of the orbit of the satellite's subsatellite point and the equator, cjIs a satellite xjThe number of the intersection periods in one day,meaning that the rounding is done down,the period of intersection of the satellite obtained in step (1) with the equatorial plane (unit is hour).
(4) And setting the sampling granularity under the full coverage condition of the equator, obtaining the longitude of the sampling point of the equator, and determining the coverage condition of the satellite load at the intersection point of the equator to the sampling point of the equator by using the satellite load breadth and all the obtained longitude of the intersection points of the equator.
Firstly, defining an equatorial sampling point longitude covering condition representation method:
carrying out rasterization equipartition sampling on the equator, wherein under the condition of N equipartition, N +1 sampling points exist, and the position of the p-th sampling point is recorded as NpTaking the coverage of the satellite width on all sampling points as a coverage condition judgment basis, and recording the number of uncovered sampling points as NIs prepared fromAnd the threshold value of the number of uncovered sampling points for finishing the calculation of the full coverage cycle is recorded as NtolIntroducing a flag function flag (p) at each sampling point,
then, determining the coverage of the satellite load at the equator intersection point on the equator sampling point:
with satellite xjFor the coverage of the c-th intersection on day d as an example, first, the satellite x is calculatedjCurrent coverage interval ofWhereinxloc(j) Is the initial position of the jth star (in radians), Δ λjIs a satellite xjAdjacent railAnd judging whether the N +1 equatorial sampling points are covered in sequence according to the longitude difference (the unit is radian) of the trace on the equator.
The step of whether the equator sampling point is covered or not is as follows: if the current sampling point NpIf the sampling point N is not covered (i.e. the flag function value flag (p) of the current sampling point is 0), the sampling point N is judged to be coveredpWhether or not in a coverage areaWithin the range: if N is presentpIn the interval, the flag function value flag (p) of the current sampling point is changed to 1, and the number of uncovered sampling points is recorded as NIs prepared from＝NIs prepared from-1; if N is presentpNot in this interval, no operation is performed. And (5) circulating for N +1 times to finish the coverage condition determination of each equatorial sampling point.
(5) And (3) taking the integral multiple of the intersection point period of the satellite in units of days, determining the initial position of the next day of the satellite and all the intersection point longitudes of the orbit of the satellite under the satellite on the day and the equator by using the steps (1) to (3), and then obtaining the coverage condition of the satellite load on the equator sampling point in one day by using the satellite load width and the equator sampling point longitude in the step (4).
The above solving process is as follows: with satellite xjTaking the coverage of c intersection points on the day d as an example, firstly, determining the coverage interval of the c intersection points by using the method for determining the coverage of the satellite load on the equator sampling point at the equator intersection point in the step (4)And (4) performing coverage judgment on the equator sampling point in the step (4) in each coverage interval in a circulating manner, so that the coverage condition of the satellite load on the equator sampling point in one day can be obtained.
(6) And (3) repeating the steps (1) to (5) by utilizing three orbit characteristic parameters of the multi-star semi-major axis, the orbit eccentricity and the orbit inclination angle to obtain the longitude of the intersection point of the equator of each star and the coverage condition of the sampling point of the equator.
The above solving process is as follows: taking the coverage situation of the day d as an example, firstly, the initial position (unit is radian) of each star and the longitude difference (unit is radian) of the adjacent track of each star on the equator, which are required by the calculation in the steps (4) and (5), are obtained by using the steps (1) to (3), and then, the step (5) is repeatedly used for each star, so that the longitude of each star equator sampling point and the coverage situation of the equator sampling point on the day d can be obtained.
(7) And (5) judging whether the equator sampling points are completely covered or not, if not, repeating the step (5), and if the equator sampling points are completely covered and the days are completely recorded, obtaining the multi-satellite combination coverage time index value.
The above solving process is as follows: assuming that the initial state d is 0, after step (6) is completed, it is determined whether to cover all the equatorial sampling points, i.e. NIs prepared fromWhether or not less than NtolIn which N isIs prepared fromNumber of uncovered samples, NtolThreshold number of uncovered samples to end the full coverage cycle, if NIs prepared fromGreater than NtolIf d is d +1, the step (6) is executed in a loop, and if N is not equal to d +1Is prepared fromIs less than or equal to NtolAnd obtaining the number of days d covered by the multi-satellite combined covering time index value.
Claims (8)
1. A system design-oriented multi-satellite combination coverage characteristic rapid determination method is characterized by comprising the following steps:
(1) determining the intersection point period of the satellite and the equatorial plane by utilizing three orbit characteristic parameters of an orbit semi-major axis, an orbit eccentricity and an orbit inclination angle;
(2) determining the longitude difference of the track of the subsatellite point of the adjacent orbit of the satellite on the equator by using the earth rotation angular velocity and the intersection point cycle obtained in the step (1);
(3) taking the intersection points of the satellite initial orbit subsatellite points and the equator as the initial positions of the satellites, and determining all the longitude of the intersection points of the satellite subsatellite point tracks and the equator in one day by the intersection point period of the satellites and the equator plane and the longitude difference of the satellite subsatellite point tracks of the adjacent orbits in the step (2) on the equator;
(4) setting sampling granularity under the full coverage condition of the equator, obtaining longitude of sampling points of the equator, and determining the coverage condition of satellite loads at the intersection points of the equator on the sampling points of the equator by using the satellite load breadth and all obtained longitude of intersection points of the equator;
(5) taking the integral multiple of the satellite intersection period and taking the number of the satellite intersection period as a unit of day, determining the initial position of the next day of the satellite and all the intersections longitude of the orbit of the satellite under the current day and the equator by using the steps (1) to (3), and obtaining the coverage condition of the satellite load on the equator sampling point in one day by using the satellite load width and the equator sampling point longitude in the step (4);
(6) repeating the steps (1) to (5) by utilizing three orbit characteristic parameters of a multi-star semi-major axis, orbit eccentricity and orbit inclination to obtain the longitude of each star equator intersection point and the coverage condition of each equator sampling point;
(7) and (5) judging whether the equator sampling points are completely covered or not, if not, repeating the step (5), and if the equator sampling points are completely covered and the days are completely recorded, obtaining the multi-satellite combination coverage time index value.
2. The method for rapidly determining the architectural design-oriented multi-satellite combined coverage property according to claim 1, wherein: the method comprises the following steps that (1) the intersection point period of the satellite and the equatorial plane is determined by utilizing three orbit characteristic parameters including orbit semimajor axis, orbit eccentricity and orbit inclination, and the intersection point period is obtained by the following method:
let the period of intersection of the satellite and the equatorial plane be
Wherein:
where n' is the satellite motion averaging period, i.e.
In the formula,is the argument of the near place;is a flat proximal angle; a is a semi-major axis of the track, and the unit is meter; e is track eccentricity, i is track inclination, J21.082625829965505e-3, second-order band harmonic coefficients of earth's gravitational potential, Re6378137.0m, the equatorial radius of the earth; mu.sE3.98600436e +14, is the gravitational constant of the earth.
3. The architectural design oriented multi-satellite combination coverage characteristic rapid determination method according to claim 2, characterized in that: the method for determining the longitude difference of the track of the subsatellite point of the adjacent orbit of the satellite in the step (2) on the equator comprises the following steps:
let satellite xjDifference in longitude Δ λ between adjacent tracks on the equatorj，
Wherein,
the rate of change of the right ascension at the intersection point;is the rotational angular velocity of the earth; a is a semi-major axis of the track; e is track eccentricity and i is track inclination.
4. The method for rapidly determining the architectural design-oriented multi-satellite combined coverage property according to claim 1, wherein: in the step (3)
(1) The method for determining the initial position of each satellite comprises the following steps:
let the initial position of the jth star be xloc(j），xloc(j) Is a random number between (0,2 pi) and is greater than or equal to the minimum load width W in the multi-star combinationminExpressed in radians; obtaining the initial position of each satellite in the same way;
(2) the longitude determination method for all intersection points of the satellite subsatellite point track and the equator in one day comprises the following steps:
Xj＝{xloc(j)+(n-1)·Δλjn belongs to N and N is less than or equal to cjIn which xloc(j) Is the initial position of the jth satellite, n is the ordinal number of the intersection of the orbit of the satellite's subsatellite point and the equator, cjIs a satellite xjThe number of the intersection periods in one day, meaning that the rounding is done down,the period of the intersection of the satellite obtained in step (1) with the equatorial plane is expressed in hours.
5. The method for rapidly determining the architectural design-oriented multi-satellite combined coverage property according to claim 1, wherein: in the step (4)
(1) The longitude coverage condition of the equatorial sampling points is expressed by the following method:
carrying out rasterization equipartition sampling on the equator, wherein under the condition of N equipartition, N +1 sampling points exist, and the position of the p-th sampling point is recorded as NpTaking the coverage of the satellite width on all sampling points as a coverage condition judgment basis, and recording the number of uncovered sampling points as NIs prepared fromAnd the threshold value of the number of uncovered sampling points for finishing the calculation of the full coverage cycle is recorded as NtolIntroducing a flag function flag (p) at each sampling point,
(2) the method for determining the coverage condition of the satellite load at the equator intersection point to the equator sampling point comprises the following steps:
for satellite xjThe coverage of the c intersection on day d is calculated by first calculating the satellite xjCurrent coverage interval ofWhereinxloc(j) Is the initial position of the jth star; delta lambdajIs a satellite xjLongitude differences of adjacent tracks on the equator are in radian units, and whether N +1 equator sampling points are covered or not is sequentially judged;
the step of judging whether the equator sampling points are covered is as follows: if the current sampling point NpIf the current sampling point is not covered, namely the flag function value flag (p) of the current sampling point is 0, the sampling point N is judgedpWhether or not in a coverage areaWithin the range: if N is presentpWithin the intervalChanging the flag function value flag (p) of the current sampling point to 1, and recording the number of uncovered sampling points as NIs prepared from＝NIs prepared from-1; if N is presentpAnd (4) not in the interval, performing no operation, circulating for N +1 times, and finishing the coverage condition determination of each equatorial sampling point.
6. The method for rapidly determining the architectural design-oriented multi-satellite combined coverage property according to claim 1, wherein: the method for determining the coverage condition of the satellite load on the equator sampling point in one day in the step (5) comprises the following steps:
for satellite xjAnd (4) determining the coverage of c intersection points on the day d by using the method for determining the coverage of the satellite load on the equator sampling point at the equator intersection point in the step (4) And (4) performing coverage judgment on the equator sampling point in the step (4) in each coverage interval in a circulating manner to obtain the coverage condition of the satellite load on the equator sampling point in one day.
7. The method for rapidly determining the architectural design-oriented multi-satellite combined coverage property according to claim 1, wherein: the method for determining the longitude and the coverage of each star equatorial sampling point in the step (6) comprises the following steps:
and (3) for the coverage condition of the day d, firstly obtaining the initial position of each satellite and the longitude difference of the adjacent track of each satellite on the equator required for calculation in the steps (4) and (5) by using the steps (1) to (3), and then repeatedly using the step (5) for each satellite, namely obtaining the longitude of each satellite equator sampling point and the coverage condition of the equator sampling point on the day d.
8. The method for rapidly determining the architectural design-oriented multi-satellite combined coverage property according to claim 1, wherein: the method for determining the number of days d covered by the multi-satellite combination coverage time index value in the step (7) comprises the following steps:
assuming that the initial state d is 0, after step (6) is completed, it is determined whether to cover all the equatorial sampling points, i.e. NIs prepared fromWhether or not less than NtolIn which N isIs prepared fromNumber of uncovered samples, NtolThreshold number of uncovered samples to end the full coverage cycle, if NIs prepared fromGreater than NtolIf d is d +1, the step (6) is executed in a loop, and if N is not equal to d +1Is prepared fromIs less than or equal to NtolAnd obtaining the number of days d covered by the multi-satellite combined covering time index value.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610972224.8A CN106570270B (en) | 2016-10-31 | 2016-10-31 | A kind of more star combined covering characteristic fast determination methods of System of System Oriented design |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610972224.8A CN106570270B (en) | 2016-10-31 | 2016-10-31 | A kind of more star combined covering characteristic fast determination methods of System of System Oriented design |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN106570270A CN106570270A (en) | 2017-04-19 |
| CN106570270B true CN106570270B (en) | 2019-09-06 |
Family
ID=58541694
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201610972224.8A Active CN106570270B (en) | 2016-10-31 | 2016-10-31 | A kind of more star combined covering characteristic fast determination methods of System of System Oriented design |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN106570270B (en) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109583055B (en) * | 2018-11-15 | 2022-02-22 | 中国人民解放军61646部队 | Satellite intersatellite point trajectory distribution optimization adjustment method based on coverage circle |
| CN111141278B (en) * | 2019-12-13 | 2021-08-10 | 航天东方红卫星有限公司 | Method for determining equatorial orbit semi-major axis by using sub-satellite point timing regression |
| CN111189458B (en) * | 2019-12-31 | 2021-04-27 | 北京跟踪与通信技术研究所 | Method and device for quickly estimating number threshold of low-orbit space accommodating satellites |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101060362A (en) * | 2007-04-13 | 2007-10-24 | 北京神州天鸿科技有限公司 | Satellite communication system-based multiple data service system |
| CN101173985A (en) * | 2006-11-01 | 2008-05-07 | 中国科学院国家天文台 | Passive radar detection method for detecting low-altitude objective by satellite signal |
| CN103678787A (en) * | 2013-11-29 | 2014-03-26 | 中国空间技术研究院 | Sub-satellite point circular geosynchronous orbit design method |
| CN105353621A (en) * | 2015-11-30 | 2016-02-24 | 北京控制工程研究所 | Fault mode thrust allocation method for geostationary orbit satellite electric thruster |
| CN106027138A (en) * | 2016-05-05 | 2016-10-12 | 清华大学 | Ground station system and method for avoiding collinear interference with geostationary satellite |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20150158602A1 (en) * | 2013-12-11 | 2015-06-11 | Tawsat Limited | Inclined orbit satellite systems |
| US9256787B2 (en) * | 2014-01-09 | 2016-02-09 | Raytheon Company | Calculation of numeric output error values for velocity aberration correction of an image |
-
2016
- 2016-10-31 CN CN201610972224.8A patent/CN106570270B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101173985A (en) * | 2006-11-01 | 2008-05-07 | 中国科学院国家天文台 | Passive radar detection method for detecting low-altitude objective by satellite signal |
| CN101060362A (en) * | 2007-04-13 | 2007-10-24 | 北京神州天鸿科技有限公司 | Satellite communication system-based multiple data service system |
| CN103678787A (en) * | 2013-11-29 | 2014-03-26 | 中国空间技术研究院 | Sub-satellite point circular geosynchronous orbit design method |
| CN105353621A (en) * | 2015-11-30 | 2016-02-24 | 北京控制工程研究所 | Fault mode thrust allocation method for geostationary orbit satellite electric thruster |
| CN106027138A (en) * | 2016-05-05 | 2016-10-12 | 清华大学 | Ground station system and method for avoiding collinear interference with geostationary satellite |
Non-Patent Citations (2)
| Title |
|---|
| 《Coverage and overlap of satellites in circular equatorial orbits. With applications to multiple-access communication-satellite systems》;S.G. Lutz;《Proceedings of the Institution of Electrical Engineers 》;19660930;第113卷(第9期);第1495-1503页 |
| 《小倾角GEO卫星多波束天线覆盖特性优化》;陈余军 等;《中国空间科学技术》;20140225(第1期);第10-17页 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN106570270A (en) | 2017-04-19 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Penny et al. | A hybrid global ocean data assimilation system at NCEP | |
| CN111045046B (en) | Short-term ionosphere forecasting method and device based on NARX | |
| CN106570270B (en) | A kind of more star combined covering characteristic fast determination methods of System of System Oriented design | |
| Elvidge et al. | Using the local ensemble transform Kalman filter for upper atmospheric modelling | |
| Gehly et al. | Sensor allocation for tracking geosynchronous space objects | |
| CN104090280B (en) | A kind of ionosphere delay correction forecasting procedure based on region CORS | |
| Han et al. | Visibility optimization of satellite constellations using a hybrid method | |
| CN109639338B (en) | Design method of global coverage constellation suitable for communication, navigation and remote integration application | |
| Habarulema | A contribution to TEC modelling over Southern Africa using GPS data | |
| Nefedyev et al. | Creation of a global selenocentric coordinate reference frame | |
| Weng et al. | Advanced data assimilation for cloud-resolving hurricane initialization and prediction | |
| Sah et al. | Constellation design of remote sensing small satellite for infrastructure monitoring in india | |
| Boniface et al. | Uncertainty quantification of the DTM2020 thermosphere model | |
| CN103235870B (en) | Take into account the sun synchronous orbit Inclination biased method of multitask height | |
| Jagoda et al. | Estimation of the Love numbers: k 2, k 3 using SLR data of the LAGEOS1, LAGEOS2, STELLA and STARLETTE satellites | |
| Renga et al. | Accurate ionospheric delay model for real-time GPS-based positioning of LEO satellites using horizontal VTEC gradient estimation | |
| Malladi et al. | Satellite constellation image acquisition problem: A case study | |
| Nag | Satellite constellation mission design using model-based systems engineering and observing system simulation experiments | |
| Xu et al. | Assimilation of high frequency radar data into a shelf sea circulation model | |
| CN106384152A (en) | PF space non-cooperative target track prediction method based on firefly group optimization | |
| Liu et al. | Guidance and control technology of spacecraft on elliptical orbit | |
| Gu et al. | A Kriging based framework for rapid satellite-to-site visibility determination | |
| CN113917505A (en) | Intelligent station surveying planning method based on orbit determination precision | |
| Chen et al. | Spacecraft autonomous GPS navigation based on polytopic linear differential inclusion | |
| Jian et al. | A GPU-based phase tracking method for planetary radio science applications |
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 |