CN106970398B - Satellite visibility analysis and ephemeris forecasting method considering satellite shielding condition - Google Patents

Satellite visibility analysis and ephemeris forecasting method considering satellite shielding condition Download PDF

Info

Publication number
CN106970398B
CN106970398B CN201710187043.9A CN201710187043A CN106970398B CN 106970398 B CN106970398 B CN 106970398B CN 201710187043 A CN201710187043 A CN 201710187043A CN 106970398 B CN106970398 B CN 106970398B
Authority
CN
China
Prior art keywords
satellite
angle
shielding
altitude
azimuth
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
CN201710187043.9A
Other languages
Chinese (zh)
Other versions
CN106970398A (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.)
PowerChina Northwest Engineering Corp Ltd
Original Assignee
PowerChina Northwest Engineering Corp Ltd
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 PowerChina Northwest Engineering Corp Ltd filed Critical PowerChina Northwest Engineering Corp Ltd
Priority to CN201710187043.9A priority Critical patent/CN106970398B/en
Publication of CN106970398A publication Critical patent/CN106970398A/en
Application granted granted Critical
Publication of CN106970398B publication Critical patent/CN106970398B/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
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/20—Integrity monitoring, fault detection or fault isolation of space segment
    • 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
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/24—Acquisition or tracking or demodulation of signals transmitted by the system
    • G01S19/27—Acquisition or tracking or demodulation of signals transmitted by the system creating, predicting or correcting ephemeris or almanac data within the receiver

Abstract

The invention discloses a satellite visibility analysis method considering satellite shielding conditions, which comprises the step 1 of measuring the shielding height of a target point satelliteAngle Ei(ii) a Step 2, resolving the satellite space position and converting coordinates; step 3, in a station center coordinate system with the ground survey station R as a coordinate origin, obtaining the altitude angle e and the azimuth angle A of the instantaneous satellite; then calculating the shielding height angle E corresponding to the time and the direction Ai(ii) a Then the satellite altitude E and the shielding altitude E of the time corresponding to the azimuth are calculatediMaking a comparison if E is less than or equal to EiIf the satellite is not visible, the satellite is rejected; if E > EiIf the satellite is visible, the satellite is reserved; the satellite is filtered according to the method, and an accurate satellite visibility result is obtained. And evaluating the accuracy factor of the single station according to the satellite visibility analysis result to realize ephemeris forecast. The method is suitable for ephemeris forecast of GNSS observation in high mountain canyon regions, urban building group effect and other regions with difficult observation, and solves the problem that the traditional ephemeris forecast result is inconsistent with the actual observation condition.

Description

Satellite visibility analysis and ephemeris forecasting method considering satellite shielding condition
Technical Field
The invention belongs to the technical field of GNSS control network measurement, and particularly relates to a satellite visibility analysis method considering a satellite shielding condition, and further relates to an ephemeris forecasting method considering the satellite shielding condition.
Background
The accuracy of the data collected is affected by the measurement accuracy of the GNSS control network, which is affected by the number of satellites that can be received in a unit time, the distribution of the satellites, the duration of received satellite signals, and other factors. Therefore, when planning a measurement scheme, ephemeris forecast is usually required to be performed on a measurement area, satellite distribution conditions and forecast accuracy factors (DOP values) in a future period of time are known in advance through the forecast ephemeris, the accuracy conditions of the whole engineering network are estimated according to the information, and then the measurement network shape and the observation scheme of the GPS are optimized to acquire more accurate measurement information.
In the ephemeris, the DOP value is determined as an important step of satellite visibility calculation, and the existing single station accuracy factor calculation satellite visibility screening generally sets a fixed altitude angle (e.g. 10 ° or 15 °) as a satellite visibility judgment condition, but actually due to the fact that the satellite position is blocked by an obstacle of a station, such as a high mountain canyon area and an urban building group effect, there is a very large difference in actually receivable satellite signals, when the cutoff altitude E is 30 °, the time for seeing 4 or more GPS satellites accounts for 90% of the whole day, and when the cutoff altitude E is 40 °, the time for seeing 4 or more GPS satellites accounts for 47% of the whole day. In high mountain canyon regions, the shielding condition of a height angle larger than 30 degrees is often encountered, which causes great errors to the satellite visibility and the accuracy of precision factor prediction.
The existing multi-satellite ephemeris forecast does not take the occlusion condition of a survey station into consideration, and the DOP value estimation result has larger deviation, so that a satellite visibility analysis method taking the occlusion condition of a satellite into consideration needs to be researched, and the accurate ephemeris forecast can be obtained only by calculating the accurate DOP value of a forecast precision factor. On the basis, the accuracy condition of the whole engineering network can be estimated more accurately, and a more optimized GNSS measurement network shape and observation scheme can be obtained.
Disclosure of Invention
The invention aims to provide a satellite visibility analysis method considering GNSS satellite shielding conditions, and solves the problem that when obstacles shield a station under a canyon region or urban building group environment, the existing GNSS satellite visibility analysis does not consider the problem that the shielding conditions of the station cause larger deviation of an analysis result.
The invention also aims to provide an ephemeris forecasting method considering the satellite shielding condition, which solves the problem of inaccurate ephemeris forecasting caused by the fact that the existing DOP value calculation does not consider the shielding condition of a measuring station when the measuring station has barrier shielding.
The invention adopts a technical scheme that a satellite visibility analysis method considering satellite shielding conditions comprises the following steps:
step 1, measuring the shielding altitude of a target point satellite;
under the station center coordinates, firstly determining the coordinates of a target point, then measuring a section at certain intervals clockwise from the north direction by taking the target point as the center, calculating the maximum shielding angle in the corresponding direction according to the measured section characteristic points, repeatedly measuring each direction, and obtaining the shielding height angle E corresponding to each direction of the target pointi。
Step 2, resolving the satellite space position and converting coordinates;
acquiring broadcast ephemeris information, then calculating the position of the satellite in the orbital plane coordinate system, and finally respectively acquiring the position of the satellite in the instantaneous spherical coordinate system and the position of the satellite in the protocol terrestrial spherical coordinate system through coordinate conversion.
Step 3, filtering the obstacle to shield the satellite
Converting the geocentric coordinates of the satellite calculated in the step 2 in the protocol terrestrial coordinate system into a station center coordinate system with the survey station as an origin, and further obtaining the altitude angle e and the azimuth angle A of the satellite; then, according to the altitude angle E and the azimuth angle A of the satellite at a certain azimuth at the designated time, and the shielding altitude angle E corresponding to each direction of the target point calculated in the step 1iObtaining the shielding height angle E of the position at the momenti(ii) a Then the satellite altitude E and the shielding altitude E of the azimuth at the moment are calculatediMaking a comparison if E is less than or equal to EiIf the satellite is not visible, the satellite is rejected; if e>EiIf the satellite is visible, the satellite is reserved; according to the method, the satellites are filtered one by one, and accurate satellite visibility results are obtained.
The technical scheme is also characterized in that:
further, the shielding height angle E in step 1iThe calculation method comprises the following steps:
the measured cross-section format is defined as: di,HiWherein D isiIs the plane distance from the feature point of the section to the target point, HiIs the elevation of the characteristic point of the section; secondly, calculating the elevation angle of the characteristic point according to the cross section
Figure BDA0001255097180000031
Setting the target point to measure n altitude angles in total, and then setting the shielding altitude angle E corresponding to each direction of the target pointiComprises the following steps:
Figure BDA0001255097180000032
further, the above
Figure BDA0001255097180000033
Is calculated by the formula
Figure BDA0001255097180000034
Or
Figure BDA0001255097180000035
In the formula H0Is the target point elevation.
Further, the method for calculating the satellite altitude e and the satellite azimuth a in step 3 comprises:
and calculating the coordinates of the satellite in a station center coordinate system with the survey station as the origin:
Figure BDA0001255097180000036
wherein the content of the first and second substances,
Figure BDA0001255097180000041
[XRYRZR]Tfor the survey station space coordinates, there are:
Figure BDA0001255097180000042
wherein, B and L are the ground longitude and latitude of the measuring station.
The altitude e corresponding to the satellite is:
Figure BDA0001255097180000043
the azimuth angle a corresponding to the satellite is:
Figure BDA0001255097180000044
in the formula:
Figure BDA0001255097180000045
further, the occlusion height angle E of the time and the orientation is specified in step 3iThe calculation method comprises the following steps:
two adjacent occlusion height angles are known: (A)i-1，Ei-1)，(Ai+1，Ei+1) And solving the shielding height angle of the designated azimuth by linear interpolation: e ═ aA + b; will (A)i-1，Ei-1)，(Ai+1，Ei+1) Substituting the above formula to obtain parameters a and b; then the satellite azimuth A of the specified time is calculatediSubstituting the formula, solving the height shielding angle in the corresponding direction as: ei＝aAi+b。
When the shielding condition is approximately linear change, the altitude angle interpolation adopts i (the shielding altitude angle quantity is set to be n, i is less than or equal to n) point linear fitting to obtain the altitude angle of the undetermined point:
E＝a0+a1A+a2A2+...+anAn
the invention adopts another technical scheme that an ephemeris forecasting method considering satellite shielding conditions comprises the following steps: the method is adopted to obtain the satellite visibility under the condition of considering the obstruction, and the direction cosine method is adopted to calculate the DOP value of the geometric precision factor based on the obtained state matrix of the observation satellite group.
The DOP value comprises a comprehensive influence precision factor GDOP, a three-dimensional position precision factor PDOP, a horizontal component precision factor HDOP, a vertical component precision factor VDOP and a clock error precision factor TDOP value.
The ephemeris forecast considering the satellite shielding condition can be obtained through the steps, and the ephemeris forecast comprises a DOP value, a satellite sky view, satellite visibility, satellite visible number and satellite elevation angle change.
The method has the advantages that the method is suitable for ephemeris forecast of GNSS observation in high mountain canyon areas, urban building group effect and other areas with difficult observation, assists operators to make an observation scheme and an observation scheduling plan, and can estimate the precision index highly consistent with the actual satellite receiving condition. The situation that the traditional ephemeris forecast index is not consistent with the reality can be solved. And (4) making an observation scheduling plan according to the forecast condition in the operation time period and the distribution condition of the measurement area. The method is suitable for GNSS observation scheduling under the high mountain canyon and urban building area effect.
Drawings
FIG. 1 is a schematic diagram of a satellite visibility analysis method and ephemeris forecast of the present invention taking into account satellite occlusion conditions;
FIG. 2 is a schematic diagram of a four-star integrated ephemeris file format.
Detailed Description
The present invention will be described in further detail with reference to the drawings and the following detailed description, but the present invention is not limited to these embodiments.
The satellite visibility analysis method considering the satellite shielding condition of the invention, referring to fig. 1, comprises the following steps:
step 1, measuring the shielding altitude of a target point satellite;
measurements of single-survey-station satellite elevation angles are made based on a topographical map, a Digital Elevation Model (DEM), or other digital geographic information products with elevation attributes.
Firstly, determining the coordinates of a target point, then measuring a section (the distance of the cross section is determined according to the actual situation) at certain intervals (for example, 5 degrees) clockwise from the north direction by taking the target point as the center, calculating the maximum shielding angle in the corresponding direction according to the measured characteristic points of the section, repeatedly measuring each direction, and obtaining the shielding height angle E corresponding to each direction of the target pointi。
The measured cross-section format is:
Di,Hi
wherein DiIs the plane distance from the feature point of the section to the target point, HiSetting target points to measure n elevation angles in total for the elevation of the section characteristic point. Elevation angle of characteristic point
Figure BDA0001255097180000061
Comprises the following steps:
Figure BDA0001255097180000062
H0for the elevation of the target point, the shielding elevation angle E corresponding to each direction of the target pointiComprises the following steps:
Figure BDA0001255097180000063
the GNSS elevation is generally expressed in the form of a standing center coordinate, and needs to be converted into the form of the standing center coordinate when the elevation measurement is performed by using a topographic map and DEM model data.
When the horizontal plane of the center-of-gravity coordinates is adopted to replace the level surface under the influence of the curvature of the earth, the calculation formula of the influence on the height difference is as follows:
Figure BDA0001255097180000064
when D is 2Km and 3Km, the effect on the height difference is 0.31m and 0.71 m. Therefore, the influence can be ignored in the general altitude measurement process, and the influence needs to be considered when higher requirements are met, and the altitude of the characteristic point
Figure BDA0001255097180000071
Becomes:
Figure BDA0001255097180000072
the height angles in all directions can be measured in sequence, and therefore the height angle measurement of the single measuring station is completed.
The altitude angle measured by a single measuring station is formatted as follows:
Figure BDA0001255097180000073
step 2, resolving the satellite space position and converting coordinates;
and (2) downloading and analyzing the latest broadcast ephemeris information according to the target point shielding altitude measured in the step (1), wherein the broadcast ephemeris downloading and analysis are the basis of the GNSS ephemeris forecasting and precision evaluating module and are the precondition for satellite position calculation.
GNSS integrated ephemeris file naming rules:
example (c): brdm3550.15p
The first four digits are the name "brdm", the next three digits are the year's day, the last digit is the time period number (typically 0), and the suffix name is the number representing the year and the file identification character "p" (p represents the GNSS four-system integrated ephemeris). The parameters and meanings of the comprehensive ephemeris are shown in table 1.
TABLE 1 ephemeris parameters and meanings
Figure BDA0001255097180000074
Figure BDA0001255097180000081
The format of the GNSS integrated ephemeris file is shown in FIG. 2.
The position of the satellite in the orbital plane coordinate system is further calculated by calculating the average angular velocity of the satellite motion, the average point angle of the satellite at the observation moment, the deviation near point angle, the true near point angle, the lift angle distance and the perturbation correction item, and finally the position of the satellite in the instantaneous spherical coordinate system and the position of the satellite in the protocol spherical coordinate system are obtained by coordinate conversion. The technical route is realized as follows:
(1) calculating the average angular velocity n of the satellite motion
Firstly, according to the parameters given in the broadcast ephemeris
Figure BDA0001255097180000082
Calculating a reference time toeAverage angular velocity n of0：
Figure BDA0001255097180000083
Wherein G is a gravitational constant, M is the total mass of the earth, and the product GM is 3.986005 × 1014m3/s2。
Then, calculating the average angular velocity n of the satellite at the observation time according to the broadcast ephemeris perturbation parameter delta n:
n＝n0+Δn
(2) calculating the mean point angle M of the satellites at the observation time
M＝M0+n(t-toe)
In the formula, M0Is a reference time toeThe average point angle in time is given by the broadcast ephemeris.
(3) Calculating the angle of approach point
The kepler equation in radians is:
E＝M+e sin E
the kepler equation expressed in degrees is:
E°＝M°+ρ°·e sin E°
the equations above can be solved by either iterative or differential correction.
(4) Calculating true proximal angle f
Figure BDA0001255097180000091
Where e is the eccentricity of the satellite orbit, given by the broadcast ephemeris. Therefore, the first and second electrodes are formed on the substrate,
Figure BDA0001255097180000092
(5) calculating the rise angle distance u'
u′＝ω+f
In the formula, ω is a near-place angular distance and is given by a broadcast ephemeris.
(6) Calculating perturbation correction term deltau、δr、δi
The following 6 perturbation parameters are given in the broadcast ephemeris: cuc，Cus，Crc，Crs，Cic，CisFrom this, the factor J can be obtained2Term-induced correction of the elevation angle uuCorrection term delta of satellite vector distance rrPerturbation correction term delta of satellite inclination angle ii. The calculation formula is as follows:
Figure BDA0001255097180000093
(7) to mu ', r', i0Performing perturbation correction
Figure BDA0001255097180000094
Wherein a is the long radius of the satellite orbit,
Figure BDA0001255097180000095
given by the broadcast ephemeris; i.e. i0Is toeThe orbit inclination angle at the moment is given by Kepler 6 parameters in the broadcast ephemeris;
Figure BDA0001255097180000101
the rate of change of i is given by the perturbation nine parameter in the broadcast ephemeris.
(8) Calculating the position of the satellite in the orbital plane coordinate system
In the orbital plane rectangular coordinate system (the origin of coordinates is located at the earth center, and the X-axis points to the elevation intersection point), the plane rectangular coordinates of the satellite are:
Figure BDA0001255097180000102
(9) calculating the longitude L of the ascending point at the observation moment
If reference is made to time toeThe right ascension at the point of time-rise intersection is
Figure BDA0001255097180000103
The rate of change of the rising point with respect to time is
Figure BDA0001255097180000104
Then the rising point right ascension Ω at the observation instant t should be:
Figure BDA0001255097180000105
Figure BDA0001255097180000106
may be given from perturbation parameters of the broadcast ephemeris.
Set the book week to openThe time of onset (0 time of Sunday) Greenwich mean sidereal time is GASTweek. Then the observation instant greenwich moment is:
GAST＝GASTweek+ωet
in the formula, ωeIs the rotational angular velocity of the earth, and has a value of omegae＝7.292115×10-5rad/s, t is the time(s) in the week, so that the longitude value of the observed instantaneous ascending intersection point can be obtained as follows:
Figure BDA0001255097180000107
order to
Figure BDA0001255097180000108
Then there are:
Figure BDA0001255097180000109
(10) computing the position of a satellite in a transient global coordinate system
Knowing the geodetic longitude L of the elevation point and the inclination i of the orbital plane, the position of the satellite in the earth-fixed coordinate system can be conveniently found by two rotations:
Figure BDA0001255097180000111
(11) calculating the position of a satellite in a protocol terrestrial coordinate system
The position of the observation instant satellite in the protocol terrestrial coordinate system is as follows:
Figure BDA0001255097180000112
and (4) resolving the satellite space position, converting coordinates, and forecasting the orbit position of the satellite in a future period of time.
And 3, filtering obstacles to shield the satellite.
And analyzing and forecasting the satellite visibility number, the visibility and the DOP value of the single survey station according to the satellite almanac and the survey station satellite altitude angle file.
And (3) converting the geocentric coordinates of the satellite in the protocol terrestrial coordinate system calculated in the step (2) into a station center coordinate system with the survey station as an origin, and further obtaining the altitude angle and the azimuth angle of the satellite, wherein the calculation methods of the altitude angle and the azimuth angle of the satellite of other systems are similar to the calculation methods. After the altitude angle and the azimuth angle of the satellite are obtained, the condition of the visible satellite at the target point position is obtained according to the shielding condition of the obstacle and the coordinate elevation information of the target point in the step 1, and the filtering of the obstacle shielding satellite and the screening of the visible satellite can be realized by combining the shielding altitude angle information of the survey station measured in the step 1. The calculation process is as follows:
(1) the calculation model of the satellite altitude and azimuth is as follows.
The coordinates of the satellite in the station center rectangular coordinate system are as follows:
Figure BDA0001255097180000113
in the formula (I), the compound is shown in the specification,
Figure BDA0001255097180000121
[XRYRZR]Tare the survey station spatial coordinates.
Figure BDA0001255097180000122
Wherein, B and L are the ground longitude and latitude of the survey station respectively.
The altitude e corresponding to the satellite is:
Figure BDA0001255097180000123
the azimuth angle a corresponding to the satellite is:
Figure BDA0001255097180000124
Figure BDA0001255097180000125
(2) screening of visible satellites
According to the altitude angle shielding information in the step 1, two adjacent shielding altitude angles are known as follows:
(Ai-1，Ei-1)，(Ai+1，Ei+1)
the occlusion height angle here is solved by linear interpolation:
E＝aA+b
will (A)i-1，Ei-1)，(Ai+1，Ei+1) The parameters a and b can be obtained by substituting the above equation. Satellite azimuth A of a given timeiSubstituting the formula, solving the height shielding angle in the corresponding direction as:
Ei＝aAi+b
comparing the altitude E and the shielding altitude E of the satelliteiJudging whether the satellite is visible or not, and judging the visible discriminant of the satellite:
e＞Ei
if E > EiIf the satellite is visible, the satellite is reserved; otherwise, the satellite is rejected.
For the shielding condition which is approximately linearly changed, the altitude angle interpolation can adopt i (the shielding altitude angle quantity is set to be n, i is less than or equal to n) point linear fitting to obtain the altitude angle of the undetermined point:
E＝a0+a1A+a2A2+...+anAn
therefore, satellite visibility is judged by filtering the conditions of the satellites blocked by the obstacles one by one.
Based on the satellite visibility analysis method considering the obstacle shielding condition, the ephemeris forecast considering the obstacle shielding condition can be obtained by further calculating the DOP value of the geometric precision factor. The method is implemented according to the following steps:
and calculating the DOP value by a direction cosine method, namely calculating by using the direction cosine of the satellite constellation, on the basis of the state matrix of the satellite group obtained by analyzing the satellite visibility method. The calculated DOP values include GDOP, PDOP, HDOP, VDOP, and TDOP values.
The technical route is as follows:
in GNSS navigation and positioning, a so-called geometric dilution of precision (DOP) value is defined as a standard for measuring the influence of the satellite space geometric distribution on the positioning precision.
The technical route of the point location observation precision evaluation submodule utilizes the coordinate of the observation station and the coordinate of the satellite to calculate the DOP value (GDOP, PDOP, HDOP, VDOP and TDOP value) between the visible satellite and the observation station to obtain the point location observation precision information.
In GNSS navigation and positioning, a so-called geometric dilution of precision (DOP) value is defined as a standard for measuring the influence of the satellite space geometric distribution on the positioning precision. When the influence of clock error precision is considered, the covariance matrix of unknown parameters is:
Figure BDA0001255097180000141
wherein each element reflects the precision information clock error precision factor TDOP under specific satellite space geometric distribution
Definition of
Figure BDA0001255097180000142
Then the corresponding error in the clock difference is
mT＝σ0·TDOP
1) Three-dimensional position accuracy factor PDOP
Definition of
Figure BDA0001255097180000143
The error in the corresponding three-dimensional position is
mP＝σ0·PDOP
2)GDOP
Combining TDOP and PDOP, a precision factor can be defined that reflects the combined effect of satellite spatial distribution on receiver clock offset and position
Figure BDA0001255097180000144
Thus, the corresponding error in the space-time accuracy is
mG＝σ0·GDOP
3) Vertical component precision factor VDOP
Definition of
Figure BDA0001255097180000145
Reflecting the influence of the space geometric distribution of the satellite on the vertical component of the position of the receiver, wherein the error in the corresponding vertical component is
mv＝σ0·VDOP
Another definition of VDOP is called elevation form factor
Figure BDA0001255097180000151
Wherein, r ═ { x, y, z } is the general location vector of the measuring station; q ═ q11,q22,q33And is a three-dimensional precision factor vector.
4) Horizontal component precision factor HDOP
Definition of
Figure BDA0001255097180000152
Through the steps, ephemeris forecast considering the satellite shielding condition can be obtained.
The method is suitable for ephemeris forecast of GNSS observation in high mountain canyon areas, urban building group effect and other areas with difficult observation, assists operators to make an observation scheme and an observation scheduling plan, and can estimate the precision index highly consistent with the actual satellite receiving condition. The problem that the traditional ephemeris forecast result is inconsistent with the actual observation condition can be solved.
The foregoing description of the invention is only a few examples, and the invention is not limited to the specific embodiments described above. The foregoing detailed description is exemplary rather than limiting in nature. All such modifications, whether made by or performed within the spirit and scope of the invention, are intended to be within the scope of the invention as defined by the appended claims.

Claims (4)

1. A satellite visibility analysis method considering satellite occlusion conditions is characterized by comprising the following steps:
step 1, measuring the shielding altitude of a target point satellite;
under the station center coordinates, firstly determining the coordinates of a target point, then measuring a section at certain intervals clockwise from the north direction by taking the target point as the center, calculating the maximum shielding angle in the corresponding direction according to the measured section characteristic points, repeatedly measuring each direction, and obtaining the shielding height angle E corresponding to each direction of the target pointi；
Step 2, resolving the satellite space position and converting coordinates;
a shielding height angle E according to step 1iAcquiring broadcast ephemeris information, then calculating the position of the satellite in a orbital plane coordinate system, and finally respectively acquiring the position of the satellite in an instantaneous spherical coordinate system and the position of the satellite in a protocol terrestrial spherical coordinate system through coordinate conversion;
step 3, filtering obstacles to shield the satellite;
converting the geocentric coordinates of the satellite calculated in the step 2 in the protocol terrestrial coordinate system into a station center coordinate system with the survey station as an origin, and further obtaining the altitude angle e and the azimuth angle A of the satellite; then, according to the altitude angle E and the azimuth angle A of the satellite at a certain azimuth at the designated time, and the shielding altitude angle E corresponding to each direction of the target point calculated in the step 1iObtaining a shielding height angle E of the position at the moment; then comparing the satellite elevation angle E of the azimuth at the moment with the shielding elevation angle E, if E is less than or equal to E, indicating that the satellite is invisible, and rejecting the satellite; if E is larger than E, the satellite is visible, and the satellite is reserved; in this way, the satellite is testedFiltering one by one to obtain accurate satellite visibility results;
the calculation method of the shielding height angle E at the specified time and direction comprises the following steps:
two adjacent occlusion height angles are known: (A)i-1，Ei-1)，(Ai+1，Ei+1) And solving the shielding height angle of the designated azimuth by linear interpolation: e ═ aA + b; will (A)i-1，Ei-1)，(Ai+1，Ei+1) Substituting the above formula to obtain parameters a and b; then the satellite azimuth A of the specified time is calculatediSubstituting the formula, solving the height shielding angle in the corresponding direction as: e ═ aAi+b；
When the shielding condition is approximately linear change, the shielding altitude angle interpolation adopts i point linear fitting to obtain the altitude angle of the undetermined point:
E＝a0+a1A+a2A2+...+anAn
in the formula, n is the number of the shielding height angles, and i is less than or equal to n.
2. The method for satellite visibility analysis considering satellite occlusion conditions according to claim 1, wherein the occlusion height angle E in step 1iThe calculation method comprises the following steps:
the measured cross-section format is defined as: di,HiWherein D isiIs the plane distance from the feature point of the section to the target point, HiIs the elevation of the characteristic point of the section; secondly, calculating the elevation angle of the characteristic point according to the cross section
Figure FDA0002410889870000021
Setting the target point to measure n altitude angles in total, and then setting the shielding altitude angle E corresponding to each direction of the target pointiComprises the following steps:
Figure FDA0002410889870000022
3. satellite visibility accounting for satellite occlusion conditions as claimed in claim 2Method of analysis, characterized in that the elevation angle of the characteristic points
Figure FDA0002410889870000023
Is calculated by the formula
Figure FDA0002410889870000024
In the formula, H0Is the target point elevation.
4. The method for analyzing the satellite visibility considering the satellite occlusion condition according to claim 1, wherein the method for calculating the satellite altitude e and the satellite azimuth A in the step 3 comprises:
acquiring coordinates of a satellite in a station center coordinate system with a survey station R as a coordinate origin:
Figure FDA0002410889870000031
in the formula (I), the compound is shown in the specification,
Figure FDA0002410889870000032
[XRYRZR]Tfor the WGS-84 coordinate vector of the station-station R, then there are:
Figure FDA0002410889870000033
wherein B and L are respectively the geodetic latitude and the geodetic longitude of the survey station R,
then the satellite altitude e is:
Figure FDA0002410889870000034
the satellite azimuth A is:
Figure FDA0002410889870000035
wherein the content of the first and second substances,
Figure FDA0002410889870000036
CN201710187043.9A 2017-03-27 2017-03-27 Satellite visibility analysis and ephemeris forecasting method considering satellite shielding condition Active CN106970398B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710187043.9A CN106970398B (en) 2017-03-27 2017-03-27 Satellite visibility analysis and ephemeris forecasting method considering satellite shielding condition

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710187043.9A CN106970398B (en) 2017-03-27 2017-03-27 Satellite visibility analysis and ephemeris forecasting method considering satellite shielding condition

Publications (2)

Publication Number Publication Date
CN106970398A CN106970398A (en) 2017-07-21
CN106970398B true CN106970398B (en) 2020-05-19

Family

ID=59329844

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710187043.9A Active CN106970398B (en) 2017-03-27 2017-03-27 Satellite visibility analysis and ephemeris forecasting method considering satellite shielding condition

Country Status (1)

Country Link
CN (1) CN106970398B (en)

Families Citing this family (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107607969B (en) * 2017-08-09 2021-01-05 东南大学 Four-system pseudo range positioning method based on DCB correction
CN109270558B (en) * 2018-09-30 2021-04-13 中国气象局气象探测中心 Shelter forecasting method for mountain base
CN109633703B (en) * 2018-12-27 2020-06-30 深圳市力合微电子股份有限公司 Beidou navigation passive positioning method for responding to sheltered scene
CN112824936A (en) * 2019-11-21 2021-05-21 百度在线网络技术（北京）有限公司 Method and device for determining height of ground object, electronic equipment and medium
CN110926474B (en) * 2019-11-28 2021-09-03 南京航空航天大学 Satellite/vision/laser combined urban canyon environment UAV positioning and navigation method
CN111505686B (en) * 2020-04-17 2021-12-31 中国科学院国家授时中心 Coarse difference elimination method based on Beidou navigation system
CN111830536A (en) * 2020-08-10 2020-10-27 南京林业大学 Satellite visibility judgment method combined with terrain influence
CN114252888A (en) * 2020-09-23 2022-03-29 苏州宝时得电动工具有限公司 Autonomous robot, base station site selection method, device and storage medium
CN112526617A (en) * 2020-11-19 2021-03-19 中国电波传播研究所（中国电子科技集团公司第二十二研究所） Ionospheric tomography system observation data simulation method based on multi-source satellite signals
CN114791299B (en) * 2021-01-25 2024-03-22 深圳市万普拉斯科技有限公司 Identification method, device and equipment for inter-layer switching of terminal equipment movement
CN113075709A (en) * 2021-03-24 2021-07-06 刘成 Vehicle-mounted satellite navigation method and device, storage medium and processor
CN115015969A (en) * 2022-06-02 2022-09-06 长安大学 GNSS satellite visibility forecasting method under mountain area sheltering environment

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9250330B2 (en) * 2007-12-07 2016-02-02 Telecommunication Systems, Inc. Method and system for selecting optimal satellites for A-GPS location of handsets in wireless networks
CN102914781A (en) * 2011-08-05 2013-02-06 北京邮电大学 Method and device for generating ephemeris message of glonass satellite signal
CN102590831B (en) * 2012-02-13 2013-11-13 北京华力创通科技股份有限公司 Method and device for judging satellite visibility under carrier rotating condition
CN103760582B (en) * 2014-01-03 2015-12-30 东南大学 A kind of optimization method blocking satellite double-difference observation structure under environment
CN104408293A (en) * 2014-11-03 2015-03-11 中国人民解放军63961部队 Satellite visibility analysis method for Beidou differential reference station layout
KR101803652B1 (en) * 2016-06-23 2017-12-29 국방과학연구소 Method and Apparatus for improving Satellite Based Augmentation System availability through geometric measuring metric development for undersampled irregularity threat model in SBAS

Also Published As

Publication number Publication date
CN106970398A (en) 2017-07-21

Similar Documents

Publication Publication Date Title
CN106970398B (en) Satellite visibility analysis and ephemeris forecasting method considering satellite shielding condition
EP3646062B1 (en) Three-dimensional city models and shadow mapping to improve altitude fixes in urban environments
CN108919634A (en) A kind of three non-non-combined observation Time Transmission system and method for difference of frequency of Beidou
CN106959456B (en) A kind of GNSS SURVEYING CONTROL NETWORK Accuracy Estimation
CN101743453A (en) The post-mission high accuracy position and azimuth determining system
CN109917494A (en) Rainfall forecast method, apparatus, equipment and storage medium
CN110907971B (en) Satellite positioning method and device for high-altitude equipment, computer equipment and storage medium
CN105044733B (en) A kind of high-precision aeronautical satellite TGD parameter calibration methods
Wright et al. The effectiveness of global positioning system electronic navigation
US20090299635A1 (en) Terrain mapping
US11913809B2 (en) Systems and methods for extending the spatial coverage of a reference pressure network
CN111830536A (en) Satellite visibility judgment method combined with terrain influence
CN116755126A (en) Beidou real-time accurate positioning method based on three-dimensional model mapping matching
CN113176596B (en) Pneumatic high-elevation constraint positioning method
CN106507915B (en) Based on high rail autonomous navigation of satellite method in weak navigation constellation signal
CN115015969A (en) GNSS satellite visibility forecasting method under mountain area sheltering environment
US20160282472A1 (en) Method and Device for Determining at least one Sample-Point-Specific Vertical Total Electronic Content
KR100448054B1 (en) Method for Preparing Geographical Information System Employing the Amended Value as Road Data
Sholarin et al. Global navigation satellite system (GNSS)
CN114088080A (en) Positioning device and method based on multi-sensor data fusion
Tran et al. Impact of the precise ephemeris on accuracy of GNSS baseline in relative positioning technique
Yakubu et al. Ramification of datum and ellipsoidal parameters on post processed differential global positioning system (DGPS) data–A case study
Dabrowski Accuracy of geopotential models used in smartphone positioning in the territory of Poland
Liu Positioning performance of single-frequency GNSS receiver using Australian regional ionospheric corrections
CN110646817A (en) Method for calculating positioning error and high-precision positioning method

Legal Events

Date Code Title Description
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