CN115015969A - GNSS satellite visibility forecasting method under mountain area sheltering environment - Google Patents
GNSS satellite visibility forecasting method under mountain area sheltering environment Download PDFInfo
- Publication number
- CN115015969A CN115015969A CN202210621360.8A CN202210621360A CN115015969A CN 115015969 A CN115015969 A CN 115015969A CN 202210621360 A CN202210621360 A CN 202210621360A CN 115015969 A CN115015969 A CN 115015969A
- Authority
- CN
- China
- Prior art keywords
- satellite
- angle
- elevation angle
- obstacle
- visibility
- 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.)
- Pending
Links
- 238000013277 forecasting method Methods 0.000 title claims abstract description 10
- 230000000007 visual effect Effects 0.000 claims description 55
- 239000011159 matrix material Substances 0.000 claims description 35
- 238000000034 method Methods 0.000 claims description 27
- 238000004364 calculation method Methods 0.000 claims description 21
- 238000012545 processing Methods 0.000 claims description 18
- 238000010586 diagram Methods 0.000 claims description 12
- 230000008569 process Effects 0.000 claims description 4
- 150000001875 compounds Chemical class 0.000 claims description 3
- 238000005259 measurement Methods 0.000 description 6
- 230000003139 buffering effect Effects 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 238000012544 monitoring process Methods 0.000 description 2
- 239000000243 solution Substances 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 239000012895 dilution Substances 0.000 description 1
- 238000010790 dilution Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
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/03—Cooperating elements; Interaction or communication between different cooperating elements or between cooperating elements and receivers
- G01S19/08—Cooperating elements; Interaction or communication between different cooperating elements or between cooperating elements and receivers providing integrity information, e.g. health of satellites or quality of ephemeris data
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
Landscapes
- Engineering & Computer Science (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Computer Security & Cryptography (AREA)
- Computer Networks & Wireless Communication (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Position Fixing By Use Of Radio Waves (AREA)
Abstract
The invention provides a GNSS satellite visibility forecasting method in a mountain area sheltering environment, which relates to the technical field of satellite visibility forecasting. The invention determines the visibility state of the satellite by judging the magnitude relation between the maximum elevation angle of the mountain obstacle and the satellite elevation angle at each azimuth angle, and circularly judges the visibility states of all satellites in the forecast time period, thereby forecasting the visibility state of the satellite under the condition of mountain area shielding.
Description
Technical Field
The invention relates to the technical field of satellite visibility forecasting, in particular to a GNSS satellite visibility forecasting method in a mountainous area sheltering environment.
Background
Since China's Beidou satellite navigation system formally provides global services in 2020, satellite constellation resources of a GNSS system are greatly enriched. When various types of GNSS measurement operations are carried out in mountainous areas, whether a GNSS receiver can receive satellite signals is crucial, and engineering construction such as mountain area search and rescue, geological disaster monitoring, roads, railways and the like all depend on the visual conditions of GNSS satellites.
At present, the visibility of satellites is mainly forecasted through satellite almanac data before GNSS measurement operation, and whether the number of visible satellites in a certain area meets the operation requirement is judged through setting a satellite altitude angle threshold (such as 5 degrees or 15 degrees). The reason why the mountain area influences the satellite visibility condition is mainly mountain body shielding, the method for judging the satellite visibility through the threshold value at the present stage cannot fully utilize the existing terrain data, the shielding condition of the mountain area terrain on the satellite is not considered, so that misjudgment on the satellite visibility state is caused, a relatively accurate satellite visibility result under the condition of the mountain area shielding cannot be obtained, and the forecasting effect is poor.
Disclosure of Invention
The invention provides a GNSS satellite visibility forecasting method in a mountain sheltered environment, and aims to overcome the defects in the prior art.
In order to achieve the above purpose, the invention provides the following technical scheme: a GNSS satellite visibility forecasting method in a mountain area sheltered environment comprises the following steps:
determining the coordinates of an operation place, a forecast time period, a forecast time interval and buffer distance parameters;
downloading Digital Elevation Model (DEM) data and satellite almanac files containing an operation site area;
performing visual field analysis on a region which takes an operation place as a center and takes a buffer distance as a radius to obtain the maximum elevation angle of the mountain obstacle at each azimuth angle of the region, and obtaining a function f of the maximum elevation angle and the azimuth angle through polynomial fitting;
the satellite visibility is judged in the traversal forecast time period, and the satellite visibility forecast result is displayed;
the process of calculating the satellite visibility through the forecast period comprises the following steps:
calculating coordinates of the satellite in a WGS-84 coordinate system and an azimuth angle and an altitude angle of the operation place and each satellite by using a satellite almanac file;
calculating the maximum elevation angle of the obstacle at each azimuth angle according to the function f by using the azimuth angle of each satellite at the current moment;
judging whether the satellite is visible at the moment according to the maximum elevation angle of the azimuth obstacle and the elevation angle of the satellite;
circularly calculating the visibility condition of each satellite in the forecast time period according to the forecast time interval;
and counting the total number of visible satellites at each moment in the forecast time period, and calculating the geometric accuracy factor of the moment.
Preferably, before obtaining the function f of the maximum elevation angle and the azimuth angle through polynomial fitting, the visibility analysis is performed on the area within the buffer distance of the operation site, and the method specifically includes the following steps:
taking the operation place as a coordinate origin, dividing the DEM into eight direction lines and 8 sectors, giving two matrixes with the same row number and column number as the DEM, respectively called as an auxiliary grid and a visual matrix, and initializing the two matrixes;
grid points on the eight direction lines are processed, and the shielding relation of the grid points is judged;
and processing the grid points in the eight sectors to complete the visual field analysis of the DEM.
Preferably, the grid point occlusion relation calculating step on the eight direction lines includes:
calculating a minimum elevation value Z of the grid point and the station which is just visible;
assigning values to the auxiliary grid and the visual matrix, if the elevation of the grid point is greater than the minimum elevation value Z, enabling the survey station and the grid point to be visual, and setting the corresponding position of the grid point in the auxiliary grid and the visual matrix as true, otherwise, enabling the grid point to be invisible and false;
and moving the grid points outwards along the direction lines, repeating the steps, and if the distance reaches the buffer distance, performing the processing of the next direction line until the eight direction lines are calculated.
Preferably, the remaining grid points that fall into the sector after the directional grid point processing are processed, and the step of processing the grid points in the sector includes:
calculating a minimum elevation value Z of the grid point and the station which is just visible;
assigning values to the auxiliary grid and the visual matrix, if the elevation of the grid point is greater than the minimum elevation value Z, enabling the survey station and the grid point to be visual, and setting the corresponding position of the grid point in the auxiliary grid and the visual matrix as true, otherwise, enabling the grid point to be invisible and false;
and moving the grid point away from the direction of the observation station, and entering the next sector if the distance reaches the buffering distance until the calculation of all eight sectors is completed, namely completing the visual field analysis of the DEM.
Preferably, extracting the maximum elevation angle of the mountain obstacle and carrying out polynomial fitting;
the method comprises the following specific steps of extracting the maximum elevation angle of the mountain obstacle:
judging any position in the visual matrix, if eight adjacent positions around the visual matrix have the invisible viewpoints, recording the position, and obtaining a position set after searching all the positions of the matrix, wherein the position set is the boundary position of the visual domain;
calculating the elevation angle of the mountain body at the corresponding position in the position set: traversing the set, and calculating the azimuth angle and the elevation angle of the position by combining the coordinates of the position in the DEM data and the coordinates of the operation place, and finally obtaining the set of the azimuth angle and the elevation angle of all boundary points, which is marked as { (A1, E1), (A2, E2), …, (An, En) };
wherein, the process of performing polynomial fitting comprises:
after extracting the maximum elevation angle of the mountain obstacle, a plurality of discrete points (A, E) can be obtained, polynomial fitting is carried out on the discrete points (A, E), the fitting criterion adopts a least square criterion, and finally, a functional relation formula E ═ f (A) of an azimuth angle and the maximum elevation angle of the obstacle is obtained, the maximum elevation angle of the mountain obstacle of the azimuth angle can be obtained through any azimuth angle, and the calculation formula comprises the following steps:
in the formula, n represents the fitting order,
representing polynomial fitting coefficients.
Preferably, the forecasting period is traversed according to a set time interval deltaT, the satellite visibility state at the moment is calculated, and the processing flow of each traversal is the same.
Preferably, the satellite visibility state calculating step is as follows:
calculating the coordinates (x) of the satellite at time t s ,y s ,z s )；
And calculating the elevation angle and the azimuth angle of each satellite sight line according to the operation place coordinates (x, y, z), wherein the specific calculation formula is as follows:
elevation＝arcsin(u)
in the above-mentioned formula, the compound of formula,
a unit vector representing a line from the operation point to the satellite view; r represents the geometric distance from the operation site to the satellite; OMGE represents the earth's angular velocity, which is a constant with a value of 7.2921151467 × 10 -5 (ii) a C is the speed of light;
a vector representing a center-of-station coordinate system centered on the operation site; azimuth, elevation represent the azimuth and elevation angles of the satellite, respectively.
Preferably, the azimuth angle of each satellite at the time t is substituted into the maximum elevation angle of the mountain obstacle and a function of the maximum elevation angle and the azimuth angle of the mountain obstacle, so that the maximum elevation angle of the obstacle of each satellite at the time t can be obtained;
comparing the altitude angle of each satellite at the time t with the calculated maximum elevation angle of the obstacle, wherein if the altitude angle of the satellite is smaller than the maximum elevation angle of the obstacle, the satellite is invisible at the time t; on the contrary, if the altitude angle of the satellite is larger than or equal to the maximum elevation angle of the obstacle, the satellite is visible at the moment t;
and after all the satellites at the time t are visually judged, counting the number of all the visible satellites at the time t.
Preferably, a geometric precision factor at time t is calculated, and the specific calculation formula includes:
wherein n is the total number of visible satellites at time t, el n El in (1) denotes the height angle, az n Az in (1) represents an azimuth, and n represents an nth satellite, e.g., el 1 Denotes the altitude, az, of the 1 st satellite 1 Representing an azimuth of the first satellite; cos and sin respectively represent cosine and sine functions; g in the G matrix represents an element thereof, wherein G 11 -g 44 Represents its position in the matrix G; GDOP, PDOP, HDOP, VDOP represent the geometric, positional, horizontal and vertical precision factors, respectively.
Preferably, the satellite visibility forecast result comprises a drawn visible satellite sky view, a satellite visibility time sequence diagram and a precision factor graph.
Compared with the prior art, the invention has the following beneficial effects:
1. the method does not need field work, directly uses the early DEM data, is convenient for programming realization, can quickly forecast the satellite visuality of a plurality of stations, does not need multiple section characteristic point measurement work, and does not increase the measurement work when the station position is changed.
2. The method for acquiring the mountain shielding maximum elevation angle after DEM visual field analysis is carried out based on the characteristics of DEM data, and different from the method based on sight shielding analysis in the scheme, the method for acquiring the mountain shielding maximum elevation angle in the invention has the advantages that the maximum elevation angle corresponding to the azimuth angle is denser, and the technical defects that the maximum shielding angle acquired by the traditional section measurement mode is sparse and the maximum shielding angle corresponding to a part of the azimuth angle is missed are overcome.
3. According to the satellite visibility state forecasting method, the satellite visibility state is determined by judging the magnitude relation between the maximum elevation angle of the mountain obstacle and the satellite elevation angle on the satellite azimuth angle, and the satellite visibility states of all satellites in a certain time period are judged in a circulating mode, so that the satellite visibility state under the condition of mountain area shielding is forecasted.
Drawings
FIG. 1 is an overall flow chart provided by the present invention;
FIG. 2 is a schematic diagram of a maximum elevation angle of a mountain obstacle according to the present invention;
FIG. 3 is a schematic diagram of the present invention providing for the grid division of DEM data;
FIG. 4 is a diagram of the result of the DEM visual domain analysis provided by the present invention;
FIG. 5 is a view of the sky of a visible satellite according to the present invention;
FIG. 6 is a timing diagram illustrating the visibility of a satellite according to the present invention;
fig. 7 is the final accuracy factor graph provided by the present invention.
Detailed Description
The following further describes embodiments of the present invention with reference to the accompanying drawings. The following examples are only for illustrating the technical solutions of the present invention more clearly, and the protection scope of the present invention is not limited thereby.
The invention provides a GNSS satellite visibility forecasting method in a mountain sheltered environment, which comprises the following steps as shown in figure 1:
step 1, determining parameters such as coordinates of a working place (a measuring station), a forecast time period, a forecast time interval, a buffer distance and the like.
And 2, downloading Digital Elevation Model (DEM) data and satellite almanac files of the operation place area.
And 3, performing visual field analysis on the operation site (the observation station) in a buffer area with a certain distance (such as 10km) to obtain the maximum elevation angle between each azimuth angle and the mountain shelter, and obtaining the functional relation between the maximum elevation angle and the azimuth angle through polynomial fitting.
And 4, calculating the coordinates of the satellite in a WGS-84 coordinate system by using the downloaded GPS, BDS, GLONASS and GALILEO satellite almanac files.
And 5, calculating the azimuth angle and the altitude angle between the operation site (the survey station) and each satellite.
And 6, calculating and obtaining the maximum elevation angle on the azimuth angle by using the functional relation obtained in the step 3, and judging whether the satellite is visible at the moment or not by comparing the size relation between the maximum elevation angle and the satellite elevation angle.
And 7, circularly calculating the visibility condition of each satellite in the forecast time period according to a certain time interval (30 s).
And 8, counting the total number of visible satellites at each moment in the forecast time period, and calculating a geometric precision factor (GDOP) at the moment.
And 9, drawing results such as a visual satellite sky view, a visual satellite number time sequence diagram, a precision factor change diagram and the like in a forecast time period.
Example 1
The method disclosed by the invention is easy to realize in programming, and compared with a single altitude angle threshold method, the forecasting result is more accurate, the method can be used for automatic site selection and other operations of GNSS disaster monitoring points in mountainous areas, and has better application value.
The "stations" appearing hereinafter are the "job sites".
Firstly, forecasting parameter determination and preparation of satellite almanac and DEM data
1. Determining forecast parameters: determining coordinates (x, y, z) of a survey station to be subjected to a satellite visual forecast, the coordinates of which are WGS-84 coordinates; giving a forecast time interval [ ts, te ], wherein the time format is (yyyy mm dd hh mm ss), and the time interval is year, month, day, hour, minute and second; determining a forecast time interval deltaT, if the given time interval is 30 seconds, calculating the satellite visibility every 30 seconds; given a buffer distance dBuff, if the given buffer distance is 10km, only the terrain in the area centered at the survey station and having a radius of 10km is considered.
2. Satellite almanac preparation: and downloading satellite almanac data issued by GPS, BDS, GLONASS and GALILEO, wherein the validity period of the satellite almanac data is supposed to contain the given forecast time period.
DEM data preparation: the DEM data can reflect the terrain of the area around the measuring station, the invention has no requirement on the source and the resolution of the used DEM data, and the invention proposes to use 30m resolution DEM data, and 30m resolution DEM data (http:// dwtkns. com/srtm30m /) provided by the American geological survey, or 12.5m resolution DEM data (https:// search. asf. alaska. edu /) provided by the Japanese space aviation research institute can be directly downloaded, wherein the measuring station needs to be positioned in the DEM data.
Second, visual field analysis of survey station and function relation of maximum elevation angle and azimuth angle of obstacle
1. Visual field analysis of the station: if a certain satellite is shielded by a mountain, the satellite is invisible; finding the maximum elevation angle of the mountain obstacle in a certain azimuth is the key to judge whether the satellite is visible, that is, if the maximum elevation angle is larger than the altitude angle of a certain satellite, the satellite is invisible, and the schematic diagram of the maximum elevation angle of the mountain obstacle is shown in fig. 2; because DEM data is grid data, the elevation angle of the sight between a mountain and the observation station is directly analyzed, interpolation processing needs to be carried out on the DEM data, and the calculation amount is large, so that the visibility in the buffer area of the observation station is analyzed by using a method without DEM interpolation calculation, and the method specifically comprises the following steps:
(1) considering the survey station (x, y, z) as the origin of coordinates, the DEM can be divided into eight direction lines and 8 sectors, as shown in fig. 3; two matrixes with the same row and column numbers as the DEM are given and respectively called as an auxiliary grid and a visual matrix, and the two matrixes are initialized: setting auxiliary grids of the survey station and eight adjacent grid points thereof as elevation values corresponding to DEM grid points, setting corresponding visual matrix point values as true, setting the other auxiliary grid point values as zero, and setting the other visual matrix point values as false, wherein at the moment, the eight adjacent grid points and the survey station are assumed to be visible; since no other grid points are shielded between the measuring station and the adjacent grid points, the assumption can be considered to be established within the accuracy range of the DEM data.
(2) And (3) processing grid points on the eight direction lines: the grid points on the eight direction lines have simple geometric characteristics, and the specific calculation step for judging the shielding relation comprises the following steps:
calculating a minimum elevation value Z of a grid point and a survey station which are just visible: for convenience of representation, the positions of the grid points and the stations in the DEM data are represented by two-dimensional coordinates, the row number and the column number of the grid points are (m, n), the row number and the column number of the stations are (i, j), and on the direction line, the point t1 between the grid points and the stations is represented as (m1, n 1); the minimum elevation Z is then calculated as:
in the above formula, subscripts denote row and column numbers in the respective grids, and r denotes an auxiliary grid; note that care should be taken to correctly specify the t1 point location; when the grid point is in the positive west direction, t1 is (m +1, n); when the grid point is in the northwest direction, t1 is (m +1, n + 1); when the grid point is located in the north direction, t1 is (m, n + 1); when the grid point is located in the northeast direction, t1 is (m-1, n + 1); and the rest directions are analogized in turn.
And (2) auxiliary grid and visual matrix assignment: if the elevation of the grid point (m, n) is greater than the minimum elevation value Z, the survey station is visible with the grid point, and the corresponding position of the grid point (m, n) in the auxiliary grid and the visual matrix is set as true; otherwise, the station is invisible to the grid point, and the corresponding position of the auxiliary grid and the visual matrix is set as false.
③ moving the lattice points (m, n) outwards along the direction line: and repeating the first step and the second step, and if the distance reaches the buffering distance, processing the next direction line until all eight direction lines are calculated.
(3) Grid points within eight sectors are processed: after the grid points on the eight direction lines are processed, the rest grid points to be processed are distributed in eight sectors, and the grid point processing step in the sectors comprises the following steps.
Calculating a minimum elevation value Z of a grid point and a survey station which are just visible: different from the grid point processing on the direction line, the grid point judgment in the sector can use two grid points t1(m1, n1) and t2(m2, n2) near the grid point as required, and the specific calculation formula comprises:
x 1 ＝m1·dx；y 1 ＝n1·dy；z 1 ＝r m1,n1 ；
x 2 ＝m2·dx；y 2 ＝n2·dy；z 2 ＝r m2,n2 ；
x 3 ＝i·dx；y 3 ＝j·dy；z 3 ＝r i,j ；
x＝m·dx；y＝n·dy；z＝Z；
x 21 ＝x 2 -x 1 ；y 21 ＝y 2 -y 1 ；z 21 ＝z 2 -z 1 ；
x 31 ＝x 3 -x 1 ；y 31 ＝y 3 -y 1 ；z 31 ＝z 3 -z 1 ；
in the above formula, r represents an auxiliary grid, subscripts of r represent row and column numbers of r, and dx and dy are the distances between the DEM grids in the direction X, Y respectively; note that care should be taken to correctly specify the t1, t2 point locations; when the grid point (m, n) is located in a sector between the west direction line and the north direction line, t1 is (m +1, n), and t2 is (m +1, n + 1); when the grid point is located in a sector between the north and the west direction lines, t1 is (m, n +1), and t2 is (m +1, n + 1); when the grid point is located in the sector between the north and northeast directions, t1 is (m-1, n +1), and t2 is (m, n + 1); and repeating the steps of the grid points in the other sectors.
And (2) auxiliary grid and visual matrix assignment: if the elevation of the grid point (m, n) is greater than the minimum elevation value Z, the survey station is visible with the grid point, and the corresponding position of the grid point (m, n) in the auxiliary grid and the visible matrix is set as true; otherwise, the station is invisible to the grid point, and the corresponding position of the auxiliary grid and the visual matrix is set as false.
Moving the grid points (m, n) away from the observation station: and (4) repeating the first step and the second step, if the distance reaches the buffer distance, entering the processing of the next sector until the calculation of all eight sectors is completed, and completing the visual field analysis of the DEM.
After the visual field analysis of the DEM data station is completed, the DEM data can be divided into 2 parts: visible and non-visible, as shown in fig. 4.
2. Extracting the maximum elevation angle of mountain shielding and fitting a polynomial:
(1) extracting the maximum elevation angle of mountain shielding: after the visual domain analysis of the DEM measuring station is carried out, a visual matrix corresponding to the size of a DEM grid can be obtained, and because the mountain visual judgment method used in the invention is used for judging based on the elevation, grid points lower than the elevation of the measuring station are also considered to be visual, namely the visual field boundary is the visible maximum range of the measuring station, the mountain elevation angle at the visual domain boundary and the mountain shielding maximum elevation angle are obtained; the specific step of extracting the maximum elevation angle of the mountain shelter comprises the following steps.
Firstly, judging any position in a visual matrix, if eight adjacent positions around the visual matrix have an invisible viewpoint, recording the position (m, n), and obtaining a position set after searching all the positions of the matrix, wherein the position set is the boundary position of the visual domain.
Calculating the elevation angle of the mountain body at the corresponding position in the position set: and traversing the set, and calculating the azimuth angle and the elevation angle of the position by combining the coordinates of the position in the DEM data and the coordinate of the station to finally obtain the set of the azimuth angle and the elevation angle of all boundary points, which is marked as { (A1, E1), (A2, E2), …, (An, En) }.
(2) Polynomial fitting: after extracting the maximum elevation angle of the mountain obstacle, a plurality of discrete points (A, E) can be obtained, polynomial fitting is carried out on the discrete points (A, E), the fitting criterion adopts a least square criterion, and finally, a functional relation formula E ═ f (A) of an azimuth angle and the maximum elevation angle of the obstacle is obtained, the maximum elevation angle of the mountain obstacle of the azimuth angle can be obtained through any azimuth angle, and the calculation formula comprises the following steps:
in the formula, n represents the fitting orderThe number of the first and second groups is counted,
representing polynomial fit coefficients.
Thirdly, calculating satellite visibility in traversal forecast period
Traversing the forecast time period [ ts, te ] according to a set time interval deltaT, calculating the satellite visibility state at the moment, wherein the processing flow of each traversal is the same, and the satellite visibility state calculation flow at a certain moment t in the forecast time period is explained below.
1. Calculating the coordinates of the satellite at the time t: the satellite almanac provides basic satellite orbit parameters and 2 clock error corrections of each satellite, and can be used for solving the approximate position and speed of each satellite, and the satellite almanac data issued by the GPS, BDS, GLONASS and GALILEO are approximately the same, taking the GPS almanac as an example, and each parameter contained therein is shown in table 1.
TABLE 1 GPS almanac parameters
The approximate position of each satellite is calculated, the calculation method by using the almanac is similar to the calculation method of the broadcast ephemeris parameters, but it is noted that the kepler orbit parameters contained in the almanac information are effective at the reference time Toe, and the ephemeris parameters not contained in the almanac can be directly zero, which has little influence on the calculation of the approximate position of the satellite; and finally, converting the obtained coordinates of each satellite under the geocentric earth-fixed state into a geodetic coordinate system WGS-84, and keeping the coordinates consistent with the coordinates of the DEM.
2. Calculating the altitude angle and azimuth angle of each satellite: obtaining WGS-84 coordinate (x) of certain satellite at time t s ,y s ,z s ) Then, the elevation angle and azimuth angle of the survey station to the satellite sight line can be calculated according to the survey station coordinates (x, y, z), and the specific calculation formula comprises:
elevation＝arcsin(u)。
in the above-mentioned formula, the compound of formula,
a unit vector representing a station-to-satellite line of sight; r represents the geometric distance from the survey station to the satellite; OMGE represents the earth's angular velocity, which is a constant with a value of 7.2921151467 × 10 -5 (ii) a C is the speed of light;
representing a station center coordinate system vector with the survey station as the center; azimuth, elevation represent the azimuth and elevation angles of the satellite, respectively.
3. Calculating the maximum elevation angle of the obstacle of each satellite at the time t: and substituting the azimuth angle of each satellite at the time t into the maximum elevation angle of the mountain obstacle and the function of the maximum elevation angle of the mountain obstacle and the azimuth angle of the mountain obstacle, so as to obtain the maximum elevation angle of the obstacle of each satellite at the time t.
4. Judging satellite visibility and counting: comparing the altitude angle of each satellite at the time t with the calculated maximum elevation angle of the obstacle, wherein if the altitude angle of the satellite is smaller than the maximum elevation angle of the obstacle, the satellite is invisible at the time t; on the contrary, if the altitude angle of the satellite is larger than or equal to the maximum elevation angle of the obstacle, the satellite is visible at the moment t; and after all the satellites at the time t are visually judged, counting the number of all the visible satellites at the time t.
5. Calculating a geometric precision factor: the Geometric dilution of precision (GDOP) is a coefficient important for measuring the positioning accuracy, and it characterizes the strength of the tracked satellite in the geometry at the time of measurement. The value of the GDOP becomes large, which causes the positioning accuracy to become poor; the small value of the GDOP can lead to the good positioning accuracy, and actually, the small GDOP means that the satellite is not concentrated in one area in space but can be uniformly distributed in different azimuth areas. In addition to the GDOP, there are accuracy factors such as PDOP, HDOP, VDOP, etc., which respectively represent a three-dimensional position accuracy factor, a planar position accuracy factor, and an elevation accuracy factor.
After the positions of all visible satellites and the positions of the stations are determined, the numerical values of all precision factors can be calculated, n is the total number of the visible satellites at the time t, and the specific calculation formula comprises the following steps:
in the above formula, el n El in (1) denotes the height angle, az n Az in (1) represents an azimuth, and n represents an nth satellite, e.g., el 1 Denotes the altitude, az, of the 1 st satellite 1 Representing an azimuth of the first satellite; cos and sin respectively represent cosine and sine functions; g in the G matrix represents an element thereof, wherein G 11 -g 44 Represents its position in the matrix G; GDOP, PDOP, HDOP, VDOP represent the geometric, positional, horizontal and vertical precision factors, respectively.
Fourthly, showing satellite visibility forecast results
After the satellite visibility state and the like are calculated at all times in the forecast period, corresponding images need to be drawn, the visibility condition of the observation station is visually displayed, and the drawing content mainly comprises the following steps:
1. visible satellite sky view: by taking the survey station as a center, the motion trail of the satellite in the forecast period can be visually displayed by utilizing the altitude angle and the azimuth angle obtained by calculating each visible satellite; wherein the range of the azimuth angle is 0-360 degrees, and the range of the elevation angle is 0-90 degrees; the bottom diagram of the sky view is composed of three circles from outside to inside and 4 straight lines, wherein the three circles represent height angles of 0 degree, 30 degrees and 60 degrees from outside to inside in sequence, and the height angle of a central point is 90 degrees; the four straight lines respectively represent the azimuth directions of true north-true south, northeast-southwest, true east-true west, and southeast-northwest. And traversing and drawing the positions of all visible satellites in the sky view at each moment in the forecast period to obtain a final visible satellite sky view, as shown in fig. 5.
2. Satellite visibility timing diagram: the satellite visibility graph can visually display the visible time period of each satellite in the forecast time period; taking each moment in the forecast period as an X axis, and taking all satellite numbers of GPS, BDS, GLONASS and GALILEO as a Y axis; a satellite visibility timing diagram can be obtained by traversing the visible satellite numbers at each time within the forecast period, as shown in fig. 6.
3. The precision factor graph is as follows: the precision factors including GDOP, PDOP, HDOP and VDOP can represent the satellite spatial configuration at a certain time, and the change condition of the satellite spatial configuration along with the time can be visually displayed by traversing the precision factor result of each time period in the forecast time period; taking each moment in a forecast period as an X axis, and taking a precision factor value as a Y axis; and traversing the precision factor results of each time in the forecast time period to obtain a final precision factor graph, as shown in fig. 7.
The above-mentioned embodiments are only preferred embodiments of the present invention, and the scope of the present invention is not limited thereto, and any simple modifications or equivalent substitutions of the technical solutions that can be obviously obtained by those skilled in the art within the technical scope of the present invention are included in the scope of the present invention.
Claims (10)
1. A GNSS satellite visibility forecasting method under a mountain area sheltered environment is characterized by comprising the following steps:
determining the coordinates of an operation place, a forecast time period, a forecast time interval and buffer distance parameters;
downloading Digital Elevation Model (DEM) data and satellite almanac files containing an operation site area;
performing visual field analysis on a region which takes an operation place as a center and takes the buffer distance as a radius to obtain the maximum elevation angle of the mountain obstacle at each azimuth angle of the region, and obtaining a function f of the maximum elevation angle and the azimuth angle through polynomial fitting;
the satellite visibility is judged in the traversal forecast time period, and the satellite visibility forecast result is displayed;
the process of calculating the satellite visibility through the forecast period comprises the following steps:
calculating the coordinates of the satellites in a WGS-84 coordinate system and the azimuth angle and the altitude angle of the operation site and each satellite by using the satellite almanac files;
calculating the maximum elevation angle of the obstacle at each azimuth angle according to the function f by using the azimuth angle of each satellite at the current moment;
judging whether the satellite is visible at the moment or not according to the magnitude relation between the maximum elevation angle of the obstacle at each azimuth angle and the elevation angle of the satellite;
circularly calculating the visibility condition of each satellite in the forecast time period according to the forecast time interval;
and counting the total number of visible satellites at each moment in the forecast time period, and calculating the geometric accuracy factor of the moment.
2. The method as claimed in claim 1, wherein before obtaining the function f of maximum elevation and azimuth through polynomial fitting, performing visual domain analysis on the region within the buffer distance of the operation site, the method comprises:
taking the operation place as a coordinate origin, dividing the DEM into eight direction lines and 8 sectors, giving two matrixes with the same row number and column number as the DEM, respectively called as an auxiliary grid and a visual matrix, and initializing the two matrixes;
grid points on the eight direction lines are processed, and the shielding relation of the grid points is judged;
and processing the grid points in the eight sectors to complete the visual field analysis of the DEM.
3. The method for GNSS satellite visibility prediction in mountain area sheltered environment according to claim 2, wherein the grid point sheltered relation calculation step in eight direction lines includes:
calculating a minimum elevation value Z just visible from the grid point and the measuring station;
assigning values to the auxiliary grid and the visual matrix, if the elevation of the grid point is greater than the minimum elevation value Z, enabling the survey station and the grid point to be visual, and setting the corresponding position of the grid point in the auxiliary grid and the visual matrix as true, otherwise, enabling the grid point to be invisible and false;
and moving the grid points outwards along the direction lines, repeating the steps, and if the distance reaches the buffer distance, performing the processing of the next direction line until the eight direction lines are calculated.
4. The GNSS satellite visibility forecasting method under the mountain sheltered environment as claimed in claim 3, characterized in that the remaining grid points falling into the sector after processing the direction grid points are processed, and the step of processing the grid points in the sector comprises:
calculating a minimum elevation value Z of the grid point and the station which is just visible;
assigning values to the auxiliary grid and the visual matrix, if the elevation of the grid point is greater than the minimum elevation value Z, enabling the survey station and the grid point to be visual, and setting the corresponding position of the grid point in the auxiliary grid and the visual matrix as true, otherwise, enabling the grid point to be invisible and false;
and moving the grid point away from the direction of the observation station, and if the distance reaches the buffer distance, entering the processing of the next sector until the calculation of all eight sectors is completed, namely completing the visual field analysis of the DEM.
5. The method according to claim 1, wherein the maximum elevation angle of mountain obstacle is extracted and polynomial fitting is performed;
the method comprises the following specific steps of extracting the maximum elevation angle of the mountain obstacle:
judging any position in the visual matrix, if eight adjacent positions around the visual matrix have the invisible viewpoints, recording the position, and after all the positions of the matrix are searched, obtaining a position set which is the boundary position of the visual field;
calculating the elevation angle of the mountain body at the corresponding position in the position set: traversing the set, and calculating the azimuth angle and the elevation angle of the position by combining the coordinates of the position in the DEM data and the coordinates of the operation site, and finally obtaining the set of the azimuth angle and the elevation angle of all boundary points, which is marked as { (A1, E1), (A2, E2), …, (An, En) };
wherein, the process of performing polynomial fitting comprises:
after extracting the maximum elevation angle of the mountain obstacle, a plurality of discrete points (A, E) can be obtained, polynomial fitting is carried out on the discrete points (A, E), the fitting criterion adopts a least square criterion, and finally the function relation formula E ═ f (A) of the azimuth angle and the maximum elevation angle of the obstacle is obtained, the maximum elevation angle of the mountain obstacle of the azimuth angle can be obtained through any azimuth angle, and the calculation formula comprises the following steps:
in the formula, n represents fittingThe order of the film is,
representing polynomial fitting coefficients.
6. The method as claimed in claim 1, wherein the prediction period is traversed according to a set time interval deltaT, and the satellite visibility state at that time is calculated, and the processing flow of each traversal is the same.
7. The method as claimed in claim 6, wherein the satellite visibility state calculating step comprises:
calculating the coordinates (x) of the satellite at time t s ,y s ,z s )；
And calculating the elevation angle and the azimuth angle of each satellite sight line according to the operation place coordinates (x, y, z), wherein the specific calculation formula is as follows:
elevation＝arcsin(u)
in the above-mentioned formula, the compound of formula,
a unit vector representing a line from the operation point to the satellite view; r represents the geometric distance from the operation site to the satellite; OMGE represents the earth's angular velocity, which is a constant with a value of 7.2921151467 × 10 -5 (ii) a C is the speed of light;
a vector representing a center-of-station coordinate system centered on the operation site; azimuth, elevation represent the azimuth and elevation angles of the satellite, respectively.
8. The method as claimed in claim 7, wherein the maximum elevation angle of the obstacle of each satellite at time t can be obtained by substituting the azimuth angle of each satellite at time t into the maximum elevation angle of the obstacle of the mountain and its function with the azimuth angle;
comparing the altitude angle of each satellite at the time t with the calculated maximum elevation angle of the obstacle, wherein if the altitude angle of the satellite is smaller than the maximum elevation angle of the obstacle, the satellite is invisible at the time t; on the contrary, if the altitude angle of the satellite is larger than or equal to the maximum elevation angle of the obstacle, the satellite is visible at the moment t;
and after all the satellites at the time t are visually judged, counting the number of all the visible satellites at the time t.
9. The method for forecasting visibility of the GNSS satellite under the sheltered environment in the mountainous area as claimed in claim 8, wherein the geometric accuracy factor at time t is calculated by the following specific calculation formula:
wherein n is the total number of visible satellites at time t, el n El in (1) denotes the height angle, az n Az in (1) represents an azimuth, and n represents an nth satellite, e.g., el 1 Denotes the altitude, az, of the 1 st satellite 1 Representing an azimuth of the first satellite; cos and sin respectively represent cosine and sine functions; g in the G matrix represents an element thereof, wherein G 11 -g 44 Represents its position in the matrix G; GDOP, PDOP, HDOP, VDOP represent the geometric, positional, horizontal and vertical precision factors, respectively.
10. The method as claimed in claim 1, wherein the satellite visibility forecast result includes a plotted view of a sky of a visible satellite, a timing diagram of a visibility of a satellite, and a precision factor graph.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210621360.8A CN115015969A (en) | 2022-06-02 | 2022-06-02 | GNSS satellite visibility forecasting method under mountain area sheltering environment |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210621360.8A CN115015969A (en) | 2022-06-02 | 2022-06-02 | GNSS satellite visibility forecasting method under mountain area sheltering environment |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN115015969A true CN115015969A (en) | 2022-09-06 |
Family
ID=83072978
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210621360.8A Pending CN115015969A (en) | 2022-06-02 | 2022-06-02 | GNSS satellite visibility forecasting method under mountain area sheltering environment |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115015969A (en) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115794414A (en) * | 2023-01-28 | 2023-03-14 | 中国人民解放军国防科技大学 | Satellite-to-ground full-view analysis method, device and equipment based on parallel computing |
| US20240118430A1 (en) * | 2022-09-27 | 2024-04-11 | Powerchina Northwest Engineering Corporation Limited | Gnss emergency monitoring error suppression method for alpine canyon complex environment |
Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP2634599A1 (en) * | 2012-02-29 | 2013-09-04 | Nxp B.V. | Satellite positioning using a sky-occlusion map |
| CN106959456A (en) * | 2017-03-27 | 2017-07-18 | 中国电建集团西北勘测设计研究院有限公司 | A kind of GNSS SURVEYING CONTROL NETWORKs Accuracy Estimation |
| CN106970398A (en) * | 2017-03-27 | 2017-07-21 | 中国电建集团西北勘测设计研究院有限公司 | Take the satellite visibility analysis and ephemeris forecasting procedure of satellite obstruction conditions into account |
| CN111275757A (en) * | 2020-01-08 | 2020-06-12 | 中国电子科技集团公司第五十四研究所 | Pseudo-satellite field simulation layout method based on DEM data processing |
| CN111596319A (en) * | 2020-05-18 | 2020-08-28 | 中国人民解放军国防科技大学 | Efficient simulation algorithm for influence of terrain occlusion on GNSS interference source action area |
| CN112162248A (en) * | 2020-08-21 | 2021-01-01 | 中国人民解放军93114部队 | Rapid calculation method and device for radar terrain shielding blind area |
| CN112558119A (en) * | 2020-11-30 | 2021-03-26 | 中航机载系统共性技术有限公司 | Satellite selection method based on self-adaptive BFO-PSO |
| US20220050211A1 (en) * | 2020-07-14 | 2022-02-17 | Spirent Communications Plc | Architecture for providing forecasts of gnss obscuration and multipath |
-
2022
- 2022-06-02 CN CN202210621360.8A patent/CN115015969A/en active Pending
Patent Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP2634599A1 (en) * | 2012-02-29 | 2013-09-04 | Nxp B.V. | Satellite positioning using a sky-occlusion map |
| CN106959456A (en) * | 2017-03-27 | 2017-07-18 | 中国电建集团西北勘测设计研究院有限公司 | A kind of GNSS SURVEYING CONTROL NETWORKs Accuracy Estimation |
| CN106970398A (en) * | 2017-03-27 | 2017-07-21 | 中国电建集团西北勘测设计研究院有限公司 | Take the satellite visibility analysis and ephemeris forecasting procedure of satellite obstruction conditions into account |
| CN111275757A (en) * | 2020-01-08 | 2020-06-12 | 中国电子科技集团公司第五十四研究所 | Pseudo-satellite field simulation layout method based on DEM data processing |
| CN111596319A (en) * | 2020-05-18 | 2020-08-28 | 中国人民解放军国防科技大学 | Efficient simulation algorithm for influence of terrain occlusion on GNSS interference source action area |
| US20220050211A1 (en) * | 2020-07-14 | 2022-02-17 | Spirent Communications Plc | Architecture for providing forecasts of gnss obscuration and multipath |
| CN112162248A (en) * | 2020-08-21 | 2021-01-01 | 中国人民解放军93114部队 | Rapid calculation method and device for radar terrain shielding blind area |
| CN112558119A (en) * | 2020-11-30 | 2021-03-26 | 中航机载系统共性技术有限公司 | Satellite selection method based on self-adaptive BFO-PSO |
Non-Patent Citations (9)
| Title |
|---|
| TIANHE REN: "《An Improved R-Index Model for Terrain Visibility Analysis for Landslide Monitoring with InSAR》", 《REMOTE SENSING》, 16 May 2022 (2022-05-16), pages 1 - 21 * |
| 吴艳兰: "《基于参考面的可视域算法》", 《测绘信息与工程》 * |
| 吴艳兰: "《基于参考面的可视域算法》", 《测绘信息与工程》, 30 March 2001 (2001-03-30), pages 19 - 21 * |
| 王可东: "《Kalman滤波基础及MATLAB仿真》", 31 January 2019, 北京航空航天大学出版社, pages: 267 - 268 * |
| 舒宝;何元浩;王利: "《一种适用于大尺度卫星导航定位基准站的网络RTK方法》", 《武汉大学学报(信息科学版)》 * |
| 舒宝;何元浩;王利: "《一种适用于大尺度卫星导航定位基准站的网络RTK方法》", 《武汉大学学报(信息科学版)》, 5 November 2021 (2021-11-05), pages 1609 - 1619 * |
| 许士杰: "《双模系统选星算法研究与全球导航性能分析》", 《中国优秀硕士学位论文全文数据库 信息科技辑》 * |
| 许士杰: "《双模系统选星算法研究与全球导航性能分析》", 《中国优秀硕士学位论文全文数据库 信息科技辑》, 15 March 2022 (2022-03-15), pages 136 - 1835 * |
| 雷虎民 等: "导弹制导控制原理（第2版）", 国防工业出版社, pages: 267 * |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20240118430A1 (en) * | 2022-09-27 | 2024-04-11 | Powerchina Northwest Engineering Corporation Limited | Gnss emergency monitoring error suppression method for alpine canyon complex environment |
| US12025716B2 (en) * | 2022-09-27 | 2024-07-02 | Powerchina Northwest Engineering Corporation Limited | GNSS emergency monitoring error suppression method for alpine canyon complex environment |
| CN115794414A (en) * | 2023-01-28 | 2023-03-14 | 中国人民解放军国防科技大学 | Satellite-to-ground full-view analysis method, device and equipment based on parallel computing |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106970398B (en) | Satellite visibility analysis and ephemeris forecasting method considering satellite shielding condition | |
| EP2149056B2 (en) | Positioning device, method and program with absolute positioning and relative positioning modes | |
| CN115015969A (en) | GNSS satellite visibility forecasting method under mountain area sheltering environment | |
| Sun et al. | Pursuing precise vehicle movement trajectory in urban residential area using multi-GNSS RTK tracking | |
| CN110261876B (en) | High-precision position-independent GNSS monitoring virtual reference method | |
| KR100496814B1 (en) | Method for obtaining road coordinates information and producing digital map using gps measurement | |
| CN110907971B (en) | Satellite positioning method and device for high-altitude equipment, computer equipment and storage medium | |
| Maciuk | GPS-only, GLONASS-only and combined GPS+ GLONASS absolute positioning under different sky view conditions | |
| CN111308457A (en) | Method, system and storage medium for north finding of pulse Doppler radar | |
| Djaja et al. | The Integration of Geography Information System (GIS) and Global Navigation Satelite System-Real Time Kinematic (GNSS-RTK) for Land use Monitoring | |
| Wright et al. | The effectiveness of global positioning system electronic navigation | |
| CN104181571A (en) | Method for rapidly measuring precision coordinate and elevation of ground point in area with weak CORS signals or without CORS signals | |
| US20090299635A1 (en) | Terrain mapping | |
| Sholarin et al. | Global navigation satellite system (GNSS) | |
| KR100448054B1 (en) | Method for Preparing Geographical Information System Employing the Amended Value as Road Data | |
| Pelc-Mieczkowska | Primary results of using hemispherical photography for advanced GPS mission planning | |
| Kamarudin et al. | Assessment on UAV onboard positioning in ground control point establishment | |
| KR100496811B1 (en) | Method for obtaining height information of buildings and producing digital map using gps measurement | |
| Davidovic et al. | Analysis of the influence of satellites constellation in gnss positioning accuracy | |
| CN117908019B (en) | Vertical line deviation resolving system and method based on radar altimeter formation measurement | |
| Kamarudin et al. | Investigation on UAV positioning inertial navigation system | |
| Mok et al. | GPS vehicle location tracking in dense high-rise environments with the minimum range error algorithm | |
| Projović | Application of relative positioning in topographic survey preparations on a full basis in artillery | |
| Vrhovski et al. | GNSS-based road user charging | |
| Opaluwa et al. | The Effect of Gps Satellite Geometry on The Precision Of DGPS Positioning İn Minna, Nigeria |
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 | ||
| RJ01 | Rejection of invention patent application after publication | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20220906 |