CN115015969A - GNSS satellite visibility forecasting method under mountain area sheltering environment - Google Patents

Abstract

The invention provides a GNSS satellite visibility forecasting method in a mountainous area sheltered environment, and relates to the technical field of satellite visibility forecasting. Domain analysis obtains the relationship between the maximum elevation angle and azimuth angle of mountain obstacles, and uses the satellite almanac file to calculate the coordinates, azimuth angle, and altitude angle of each satellite. Whether the moment is visible or not, the visibility forecast is carried out within the forecast time period. The invention determines the visibility state of the satellite by judging the relationship between the maximum elevation angle of the mountain obstacle on each azimuth angle and the altitude angle of the satellite, and cyclically judges the visibility states of all satellites in the forecast period, thereby forecasting the occlusion situation in the mountainous area Compared with the single altitude threshold method, the forecast results are more accurate.

Description

A GNSS satellite visibility forecasting method in mountainous area sheltered environment

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 sheltered environment.

Background technique

Since my country's Beidou satellite navigation system officially provided global services in 2020, it has greatly enriched the satellite constellation resources of the GNSS system. When carrying out various types of GNSS surveying operations in mountainous areas, it is very important whether the GNSS receiver can receive satellite signals, such as mountain search and rescue, geological disaster monitoring, road, railway and other engineering construction, all rely on the visibility of GNSS satellites.

At present, satellite visibility is mainly predicted by satellite almanac data before GNSS measurement operations, and whether the number of visible satellites in a certain area can meet the operational needs is determined by setting a satellite altitude angle threshold (such as 5° or 15°). The main reason why mountainous areas affect the visibility of satellites is the occlusion of mountains. At this stage, the method of judging satellite visibility through thresholds fails to make full use of the existing terrain data, and does not consider the occlusion of satellites by mountainous terrain, resulting in the visibility of some satellites. It is impossible to obtain more accurate satellite visibility results in the case of mountain occlusion, and the forecast effect is poor.

SUMMARY OF THE INVENTION

The present invention provides a GNSS satellite visibility forecasting method in a mountainous area sheltered environment, and aims to solve the shortcomings existing in the prior art.

In order to achieve the above object, the present invention provides the following technical solutions: a GNSS satellite visibility forecasting method in a mountainous area occluded environment, comprising the following steps:

Determine the coordinates of the operation site, forecast time period, forecast time interval and buffer distance parameters;

Download digital elevation model DEM data and satellite almanac files containing the operating site area;

Perform visual field analysis on the area with the operating site as the center and the buffer distance as the radius, obtain the maximum elevation angle of the mountain obstacles at each azimuth angle in the area, and obtain the function f of the maximum elevation angle and azimuth angle through polynomial fitting;

Traverse the forecast period to judge the satellite visibility, and display the satellite visibility forecast results;

The process of traversing the forecast period to calculate the satellite visibility includes the following steps:

Use the satellite almanac file to calculate the coordinates of the satellite in the WGS-84 coordinate system and the azimuth and altitude of the operating site and each satellite;

Using the azimuth angle of each satellite at the current moment, calculate the maximum elevation angle of each azimuth angle obstacle according to the function f;

Determine whether the satellite is visible at this moment by the size of the maximum elevation angle of the azimuth obstacle and the satellite elevation angle;

Circularly calculate the visibility status of each satellite within the forecast time period according to the forecast time interval;

Count the total number of visible satellites at each moment in the forecast period, and calculate the geometric precision factor at that moment.

Preferably, before obtaining the function f of the maximum elevation angle and azimuth angle by polynomial fitting, the visibility analysis is performed on the area within the buffer distance of the work site, and the specific steps include:

Considering the work site as the origin of coordinates, the DEM can be divided into eight direction lines and eight sectors. Given two matrices with the same row and column numbers as the DEM, they are called auxiliary grids and visual matrices, respectively. The matrix is initialized;

Process the grid points on the eight direction lines and judge their occlusion relationship;

The grid points within the eight sectors are processed to complete the visual field analysis of the DEM.

Preferably, the step of calculating the occlusion relationship of grid points on the eight direction lines includes:

Calculate the minimum elevation value Z that is just visible to the grid point and the station;

Assign values to the auxiliary grid and the visible matrix. If the grid point elevation is greater than the minimum elevation value Z, the station and grid point are visible, and the corresponding position of the grid point in the auxiliary grid and visible matrix is set to true, otherwise it is not possible. sight and false;

Move the grid point outward along the direction line, repeat the above steps, if the distance reaches the buffer distance, enter the processing of the next direction line, until the calculation of the eight direction lines is completed.

Preferably, after the direction line grid points are processed, the remaining grid points falling in the sector are processed, and the processing steps of the grid points in the sector include:

Calculate the minimum elevation value Z that is just visible to the grid point and the station;

Assign values to the auxiliary grid and the visible matrix. If the grid point elevation is greater than the minimum elevation value Z, the station and grid point are visible, and the corresponding position of the grid point in the auxiliary grid and visible matrix is set to true, otherwise it is not possible. sight and false;

Move the grid point away from the station. If the distance reaches the buffer distance, the next sector will be processed until the calculation of all eight sectors is completed, that is, the DEM visual field analysis is completed.

Preferably, extract and perform polynomial fitting on the maximum elevation angle of mountain obstacles;

Among them, the specific steps of extracting the maximum elevation angle of the mountain obstacle include:

Judging any position in the visible matrix, if there are invisible points in eight adjacent positions around it, record the position. After searching all positions in the matrix, a position set is obtained, which is the boundary position of the visual field. ;

Calculate the elevation angle of the mountain at the corresponding position in the position set: traverse the set, calculate the azimuth and elevation angles of the coordinates of the position in the DEM data and the coordinates of the work site jointly, and finally obtain the azimuth and elevation angles of all boundary points, denoted as { (A1,E1),(A2,E2),…,(An,En)};

Among them, the process of performing polynomial fitting includes:

After extracting the maximum elevation angle of the obstacles on the mountain, multiple discrete points (A, E) can be obtained, and polynomial fitting is performed on them. The fitting criterion adopts the least squares criterion, and finally the functional relationship between the azimuth angle and the maximum elevation angle of the obstacle is obtained E = f(A), through any azimuth angle, the maximum elevation angle of the mountain obstacle at this azimuth angle can be obtained. The calculation formula includes:

表示多项式拟合系数。In the formula, n represents the fitting order,

Represents the polynomial fit coefficients.

Preferably, the forecast period is traversed according to the set time interval deltaT, and the satellite visibility state at this moment is calculated, and the processing flow of each traversal is the same.

Preferably, the satellite visibility state calculation steps are as follows:

Calculate the coordinates of the satellite at time t (x ^s , y ^s , z ^s );

Calculate the altitude angle and azimuth angle of each satellite line of sight according to the coordinates of the operation site (x, y, z). The specific calculation formula is as follows:

elevation=arcsin(u)

表示作业地点至卫星视线向的单位向量；r表示作业地点至卫星的几何距离；OMGE表示地球角速度，为一常数，其值为7.2921151467×10^-5；C为光速；

表示以作业地点为中心的站心坐标系向量；azimuth、elevation分别表示卫星的方位角与高度角。In the above formula,

Represents the unit vector from the operating site to the satellite line of sight; r represents the geometric distance from the operating site to the satellite; OMGE represents the angular velocity of the earth, which is a constant, and its value is 7.2921151467×10 ^-5 ; C is the speed of light;

Represents the vector of the station center coordinate system centered on the operating site; azimuth and elevation represent the azimuth and elevation angles of the satellites, respectively.

Preferably, the azimuth angle of each satellite at time t is substituted into the maximum elevation angle of the mountain obstacle and its function with the azimuth angle, and the maximum elevation angle of the obstacle at time t of each satellite can be obtained;

Compare the altitude angle of each satellite at time t with the calculated maximum elevation angle of the obstacle. If the altitude angle of the satellite is less than the maximum elevation angle of the obstacle, the satellite is invisible at time t; on the contrary, if the altitude angle of the satellite is greater than or equal to the maximum elevation angle of the obstacle, then This satellite is visible at time t;

After the visibility of all satellites at time t is judged, the number of all visible satellites at time t is counted.

Preferably, the geometric precision factor at time t is calculated, and the specific calculation formula includes:

Among them, n is the total number of visible satellites at time t, el in el _n represents the altitude angle, az in az _n represents the azimuth angle, n represents the nth satellite, such as el ₁ represents the altitude angle of the first satellite, az ₁ represents the azimuth of the first satellite; cos and sin represent cosine and sine functions respectively; g in the G matrix represents its element, where g ₁₁ -g ₄₄ represents its position in matrix G; GDOP, PDOP, HDOP, VDOP represents geometric precision factor, position precision factor, horizontal precision factor and vertical precision factor, respectively.

Preferably, the satellite visibility prediction result includes a drawn visible satellite sky view, a satellite visibility time sequence diagram and a precision factor diagram.

Compared with the prior art, the present invention has the following beneficial effects:

1. The present invention does not require field work, and directly uses the DEM data in the early stage, which is convenient for programming, can quickly predict the satellite visibility of multiple stations, does not require multiple cross-section feature point measurement work, and does not increase when changing the position of the station. measurement work.

2. The method of the present invention to obtain the maximum elevation angle of mountain occlusion after the DEM visual field analysis is carried out based on the characteristics of the DEM data. Different from the method based on the line of sight occlusion analysis of the above scheme, the azimuth angle obtained by the method of the present invention corresponds to the maximum elevation angle. It is more dense, which solves the technical defect that the maximum occlusion angle obtained by the traditional cross-sectional measurement method is sparse, and the maximum occlusion angle corresponding to a part of the azimuth angle is missed.

3. The present invention determines the visibility state of the satellite by judging the relationship between the maximum elevation angle of the mountain obstacle on the satellite azimuth and the satellite elevation angle, and cyclically judges the visibility states of all satellites within a certain period of time, thereby forecasting. Compared with the single altitude angle threshold method, the forecast results of the satellite visibility state under the occlusion of mountainous areas are more accurate and have better application value.

Description of drawings

Fig. 1 is the overall flow chart provided by the present invention;

2 is a schematic diagram of the maximum elevation angle of a mountain obstacle provided by the present invention;

3 is a schematic diagram of dividing a DEM data grid provided by the present invention;

Fig. 4 is a DEM visual field analysis result diagram provided by the present invention;

Fig. 5 is the visible satellite sky view provided by the present invention;

Fig. 6 is the satellite visibility timing chart provided by the present invention;

FIG. 7 is the final precision factor diagram provided by the present invention.

Detailed ways

The specific embodiments of the present invention will be further described below with reference to the accompanying drawings. The following examples are only used to illustrate the technical solutions of the present invention more clearly, and cannot be used to limit the protection scope of the present invention.

The present invention provides a GNSS satellite visibility prediction method in a mountainous area sheltered environment, as shown in FIG. 1 , comprising the following steps:

Step 1: Determine the coordinates of the operation site (station), forecast time period, forecast time interval, buffer distance and other parameters.

Step 2: Download digital elevation model (DEM) data and satellite almanac files including the operating site area.

Step 3: Perform visual field analysis on the operation site (station) within a certain distance (such as 10km) buffer area to obtain the maximum elevation angle between the azimuth angle and the mountain occluder, and obtain the maximum elevation angle and azimuth angle through polynomial fitting. functional relationship.

Step 4, using the downloaded GPS, BDS, GLONASS, GALILEO satellite almanac files, calculate the coordinates of the satellites in the WGS-84 coordinate system.

Step 5: Calculate the azimuth angle and elevation angle of the operation site (station) and each satellite.

Step 6, use the functional relationship obtained in step 3 to obtain the maximum elevation angle on the azimuth angle, and judge whether the satellite is visible at this moment by comparing the magnitude relationship between the maximum elevation angle and the satellite elevation angle.

Step 7: Circularly calculate the visibility status of each satellite in the forecast time period at a certain time interval (30s).

Step 8: Count the total number of visible satellites at each moment in the forecast time period, and calculate the geometrical precision factor (GDOP) at the moment.

Step 9: Draw the results of the visible satellite sky view, the time sequence diagram of the number of visible satellites, and the variation diagram of the precision factor within the forecast period.

Example 1

The present invention uses the GNSS satellite almanac data to calculate the approximate position of the satellite, calculates the azimuth and altitude of the line connecting a certain operating site and the satellite in the mountainous area, and uses the digital elevation model information (Digital Elevation Model, DEM) of the operating mountainous area to simulate the actual terrain of the area , by judging the relationship between the maximum elevation angle of the mountain obstacle on the satellite azimuth and the satellite altitude angle, to determine the visibility state of the satellite, and cyclically determine the visibility state of all satellites in a certain period of time, so as to predict the occlusion of the mountainous area. The method disclosed in this patent is easy to program and implement, and compared with the single altitude angle threshold method, the prediction results are more accurate, and it can be used for automatic site selection of GNSS monitoring points for geological hazards in mountainous areas. good application value.

The "station" that appears below is the "operating site".

1. Determination of forecast parameters and preparation of satellite almanac and DEM data

1. Determination of forecast parameters: Determine the coordinates (x, y, z) of the station to be used for satellite visual forecast, and the coordinate system of this coordinate is the WGS-84 coordinate system; given the forecast time interval [ts, te], the time The format is (yyyy mm dd hh mm ss), which are year, month, day, hour, minute, and second respectively; determine the forecast time interval deltaT, if the given time interval is 30 seconds, the satellite visibility is calculated every 30 seconds ; Given the buffer distance dBuff, if the given buffer distance is 10km, only the terrain in the area with the station as the center and the radius of 10km is considered.

2. Satellite almanac preparation: Download the satellite almanac data released by GPS, BDS, GLONASS, and GALILEO. The validity period of the satellite almanac data should include the given forecast time period.

3. DEM data preparation: DEM data can reflect the topography of the surrounding area of the station. The present invention basically does not require the source and resolution of the DEM data used. It is recommended to use 30m resolution DEM data, which can be directly downloaded from the US Geological Survey. 30m resolution DEM data (http://dwtkns.com/srtm30m/), or 12.5m resolution DEM data provided by JAXA (https://search.asf.alaska.edu/), where the station Need to be in DEM data.

2. Analysis of the visual field of the station and the functional relationship between the maximum elevation angle and azimuth angle of obstacles

1. Analysis of the visual field of the station: if a satellite is blocked by a mountain, the satellite cannot be seen; finding the maximum elevation angle of the mountain obstacle at a certain azimuth is the key to judging whether the satellite is visible, that is, if the maximum elevation angle is greater than a certain satellite. If the altitude angle is high, the satellite is not visible, and the schematic diagram of the maximum elevation angle of the mountain obstacle is shown in Figure 2; since the DEM data is grid data, it is necessary to directly analyze the elevation angle of the line of sight between the mountain and the station, and the DEM data needs to be interpolated. The amount of calculation is large, so this patent uses a method that does not require DEM interpolation calculation to analyze the visibility in the buffer area of the station. The specific steps include:

(1) Considering the station (x, y, z) as the coordinate origin, the DEM can be divided into eight direction lines and eight sectors, as shown in Figure 3; given two matrices with the same row and column numbers as the DEM , respectively called auxiliary grid and visual matrix, and initialize these two matrices: set the auxiliary grid of the station and its eight adjacent grid points as the elevation values of the corresponding DEM grid points, and the corresponding visual The value of matrix points is set to true, the values of other auxiliary grid points are zero, and the values of other visible matrix points are set to false. At this time, it is assumed that eight adjacent grid points are visible to the station; There is no other grid point occlusion between them, and this assumption can be considered to be true within the accuracy range of DEM data.

(2) Handling the grid points on the eight direction lines: the geometric characteristics of the grid points located on the eight direction lines are simple, and the specific calculation steps for judging the occlusion relationship include:

①Calculate the minimum elevation value Z that is just visible to the grid point and the station: For convenience, two-dimensional coordinates are used to represent the position of the grid point and the station in the DEM data. The row and column numbers are (i, j), and on the direction line, the point t1 between the grid point and the station is represented as (m1, n1); the calculation formula of the minimum elevation Z is:

In the above formula, the subscript represents the row and column number in the corresponding grid, and r represents the auxiliary grid; it should be noted that the position of the t1 point should be specified correctly; when the grid point is in the due 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 in the north direction, t1 is (m, n+1); when the grid point is in the northeast direction, t1 is (m -1,n+1); other directions can be deduced and so on.

②Assignment of auxiliary grid and visual matrix: If the elevation of grid point (m,n) is greater than the minimum elevation value Z, the station and grid point are visible, and the grid point (m,n) is placed in the auxiliary grid and visual matrix The corresponding position in the grid is set to true; otherwise, the station and grid points are not visible, and the corresponding position in the auxiliary grid and visible matrix is set to false.

③ Move the grid point (m, n) outward along the direction line: repeat steps ① and ②, if the distance reaches the buffer distance, enter the processing of the next direction line, until the calculation of the eight direction lines is completed.

(3) Processing the grid points in the eight sectors: After processing the grid points on the eight direction lines, the remaining grid points to be processed are distributed in the eight sectors, and the processing steps for the grid points in the sectors include.

①Calculate the minimum elevation value Z that is just visible between the grid point and the station: Different from the grid point processing on the direction line, the grid point judgment in the sector can use the two grid points t1 (m1, n1) near the grid point as needed. and t2(m2,n2), the specific calculation formula includes:

x ₁ =m1·dx; y ₁ =n1·dy; z ₁ =r _m1,n1 ;

x ₂ =m2·dx; y ₂ =n2·dy; z ₂ =r _m2,n2 ;

x ₃ =i·dx; y ₃ =j·dy; z ₃ =r _i,j ;

x=m·dx; y=n·dy; z=Z;

x ₂₁ =x ₂ -x ₁ ; y ₂₁ =y ₂ -y ₁ ; z ₂₁ =z ₂ -z ₁ ;

x ₃₁ =x ₃ -x ₁ ; y ₃₁ =y ₃ -y ₁ ; z ₃₁ =z ₃ -z ₁ ;

In the above formula, r represents the auxiliary grid, the subscript of r represents its row and column number, and dx and dy are the spacing of the DEM grid in the X and Y directions respectively; it should be noted that the positions of t1 and t2 should be correctly specified. ;When the grid point (m,n) is located in the sector between the west and northwest direction lines, t1 is (m+1,n), t2 is (m+1,n+1); when the grid point is located in the northwest to north When the sector between the direction lines, t1 is (m, n+1), t2 is (m+1, n+1); when the grid point is located in the sector between the north and northeast direction lines, t1 is (m -1, n+1), t2 is (m, n+1); the grid points in the other sectors can be deduced in turn.

②Assignment of auxiliary grid and visual matrix: If the elevation of grid point (m,n) is greater than the minimum elevation value Z, the station and grid point are visible, and the grid point (m,n) is placed in the auxiliary grid and visual matrix The corresponding position in the grid is set to true; otherwise, the station and grid points are not visible, and the corresponding position in the auxiliary grid and visible matrix is set to false.

③ Move the grid point (m, n) away from the station: repeat steps ① and ②, if the distance reaches the buffer distance, enter the processing of the next sector until the calculation of all eight sectors is completed, that is, the visualization of the DEM is completed. Domain Analysis.

After completing the visibility area analysis of the DEM data station, the DEM data can be divided into two parts: the visible area and the invisible area, as shown in Figure 4.

2. Mountain occlusion maximum elevation angle extraction and polynomial fitting:

(1) Extraction of the maximum elevation angle occluded by the mountain: After analyzing the visual field of the DEM station, a visual matrix corresponding to the size of the DEM grid can be obtained. Since the visual judgment method of the mountain used in this patent is based on the elevation, The grid points lower than the elevation of the station are also considered visible, that is, the boundary of the viewing area is the maximum range that the station can see, so the elevation angle of the mountain at the boundary of the visible area and the maximum elevation angle of the mountain occlusion; extract the maximum elevation angle of the mountain occlusion. Specific steps include.

① Judging any position in the visible matrix, if there are invisible points in eight adjacent positions around it, record the position (m, n), and after searching all the positions of the matrix, a position set is obtained, and the set is The location of the boundaries of the viewport.

②Calculate the elevation angle of the mountain at the corresponding position in the position set: traverse the set, calculate the azimuth and elevation angles of the coordinates of the position in the DEM data and the coordinates of the station jointly, and finally obtain the azimuth and elevation angles of all boundary points, denoted as {(A1,E1),(A2,E2),…,(An,En)}.

(2) Polynomial fitting: After extracting the maximum elevation angle of the obstacles on the mountain, multiple discrete points (A, E) can be obtained, and polynomial fitting is performed on them. The functional relationship of the elevation angle E=f(A), through any azimuth angle, the maximum elevation angle of the mountain obstacle at the azimuth angle can be obtained. The calculation formula includes:

表示多项式拟合系数。In the formula, n represents the fitting order,

Represents the polynomial fit coefficients.

3. Traverse the forecast period to calculate the satellite visibility

The forecast period [ts,te] is traversed according to the set time interval deltaT, and the satellite visibility status at this moment is calculated. The processing flow of each traversal is the same. The calculation process of sexual state is explained.

1. Calculate the coordinates of the satellites at time t: The satellite almanac provides the basic satellite orbit parameters and 2 clock error correction numbers of each satellite, which can be used to solve the approximate position and speed of each satellite, published by GPS, BDS, GLONASS, GALILEO The satellite almanac data is roughly the same, taking the GPS almanac as an example, the parameters contained in it are shown in Table 1.

Table 1 GPS almanac parameters

To calculate the approximate position of each satellite, the almanac calculation is similar to the broadcast ephemeris parameter calculation method, but it should be noted that the Kepler orbit parameters contained in the almanac information are valid at the reference time Toe, and the ephemeris not included in the almanac information The parameter can be directly set to zero, which has little effect on the calculation of the approximate position of the satellite; finally, the obtained coordinates of each satellite at the center of the earth are converted into the geodetic coordinate system WGS-84, which is consistent with the coordinates of the DEM.

2. Calculate the altitude and azimuth of each satellite: After obtaining the WGS-84 coordinates (x ^s , y ^s , z ^s ) of a satellite at time t, the station can be calculated according to the station coordinates (x, y, z) The altitude angle and azimuth angle to the satellite line of sight, the specific calculation formula includes:

elevation=arcsin(u).

表示测站至卫星视线向的单位向量；r表示测站至卫星的几何距离；OMGE表示地球角速度，为一常数，其值为7.2921151467×10^-5；C为光速；

表示以测站为中心的站心坐标系向量；azimuth、elevation分别表示卫星的方位角与高度角。In the above formula,

Represents the unit vector of the line-of-sight direction from the station to the satellite; r represents the geometric distance from the station to the satellite; OMGE represents the angular velocity of the earth, which is a constant, and its value is 7.2921151467×10 ^-5 ; C is the speed of light;

Represents the station center coordinate system vector centered on the station; azimuth and elevation represent the azimuth angle and elevation angle of the satellite, respectively.

3. Calculate the maximum elevation angle of the obstacles of each satellite at time t: Substitute the azimuth angle of each satellite at time t into the maximum elevation angle of the mountain obstacle and its function with the azimuth angle, and the maximum elevation angle of each satellite at time t can be obtained. .

4. Judge the satellite visibility and make statistics: Compare the altitude angle of each satellite at time t with the calculated maximum elevation angle of the obstacle. If the satellite altitude angle is less than the maximum elevation angle of the obstacle, the satellite will not be visible at time t; otherwise, if the satellite is not visible at time t If the altitude angle is greater than or equal to the maximum elevation angle of the obstacle, the satellite is visible at time t; after judging the visibility of all satellites at time t, the number of all visible satellites at time t is counted.

5. Calculate the geometric precision factor: Geometric dilution of precision (GDOP) is a very important coefficient to measure the positioning accuracy, which characterizes the strength of the tracked satellite in the geometric structure during the measurement. The larger the value of GDOP, the worse the positioning accuracy; the smaller the value of GDOP, the better the positioning accuracy. In fact, a smaller GDOP means that the satellites are not concentrated in one area in space, but can be distributed in different areas. The azimuth area is evenly distributed. In addition to GDOP, there are PDOP, HDOP, VDOP and other precision factors, which respectively represent the three-dimensional position precision factor, the plane position precision factor and the elevation precision factor.

After determining the positions of all visible satellites and stations, the value of each precision factor can be calculated, and n is the total number of visible satellites at time t. The specific calculation formula includes:

In the above formula, el in el _n represents the elevation angle, az in az _n represents the azimuth angle, and n represents the nth satellite, for example, el ₁ represents the elevation angle of the first satellite, and az ₁ represents the azimuth of the first satellite. Angle; cos and sin represent cosine and sine functions respectively; g in G matrix represents its element, where g ₁₁ -g ₄₄ represents its position in matrix G; GDOP, PDOP, HDOP, VDOP represent geometric precision factor, position Precision factor, horizontal precision factor, and vertical precision factor.

4. Display of satellite visibility forecast results

After completing the calculation of the satellite visibility status at all times in the forecast period, it is necessary to draw corresponding images to visually display the visibility of the station. The drawing contents mainly include:

1. Visible satellite sky view: With the station as the center, the altitude and azimuth angles obtained by calculating each visible satellite can be used to intuitively display the trajectory of the satellite during the forecast period; the azimuth angle ranges from 0 to 360°, The altitude angle ranges from 0 to 90°; the base image of the sky view is composed of three circles and four straight lines from the outside to the inside. The altitude angle is 90°; the four straight lines represent the azimuth directions of due north-due south, northeast-southwest, due east-due west, and southeast-northwest respectively. By traversing and drawing the positions of all visible satellites in the sky view at each moment in the forecast period, the final visible satellite sky view is obtained, as shown in Figure 5.

2. Satellite visibility time sequence diagram: The satellite visibility diagram can intuitively display the visibility period of each satellite in the forecast period; it takes each moment in the forecast period as the X-axis, and all satellites of GPS, BDS, GLONASS and GALILEO The number is the Y-axis; by traversing the visible satellite numbers at each moment in the forecast period, the satellite visibility time sequence diagram can be obtained, as shown in Figure 6.

3. Precision factor graph: The precision factor includes GDOP, PDOP, HDOP, and VDOP, which can characterize the quality of the satellite space configuration at a certain time. By traversing the precision factor results of each time period in the forecast period, it can intuitively display the satellite space configuration over time. It takes each moment in the forecast period as the X-axis, and the precision factor value as the Y-axis; traversing the results of various precision factors at each moment in the forecast period, the final precision factor map can be obtained, as shown in Figure 7.

The above-mentioned embodiments are only preferred specific embodiments of the present invention, and the protection scope of the present invention is not limited thereto. Any person skilled in the art can obviously obtain the simplicity of the technical solution within the technical scope disclosed in the present invention. Changes or equivalent replacements all belong to the protection scope of the present invention.

Claims (10)

1. a GNSS satellite visibility forecasting method under a mountainous area sheltering environment, is characterized in that, comprises the steps: Determine the coordinates of the operation site, forecast time period, forecast time interval and buffer distance parameters; Download digital elevation model DEM data and satellite almanac files containing the operating site area; Perform visual field analysis on the area with the operating site as the center and the buffer distance as the radius, obtain the maximum elevation angle of the mountain obstacles at each azimuth angle in the area, and obtain the function f of the maximum elevation angle and azimuth angle through polynomial fitting; Traverse the forecast period to judge the satellite visibility, and display the satellite visibility forecast results; The process of traversing the forecast period to calculate the satellite visibility includes the following steps: Use the satellite almanac file to calculate the coordinates of the satellite in the WGS-84 coordinate system and the azimuth and altitude of the operating site and each satellite; Using the azimuth angle of each satellite at the current moment, calculate the maximum elevation angle of each azimuth angle obstacle according to the function f; Judging whether the satellite is visible at this moment by the relationship between the maximum elevation angle of the obstacle at each azimuth angle and the altitude angle of the satellite; Circularly calculate the visibility status of each satellite within the forecast time period according to the forecast time interval; Count the total number of visible satellites at each moment in the forecast period, and calculate the geometric precision factor at that moment. 2. the GNSS satellite visibility forecasting method under a kind of mountainous area shielding environment as claimed in claim 1, is characterized in that, before obtaining the function f of maximum elevation angle and azimuth angle by polynomial fitting, to the operation site buffering distance. The specific steps include: Considering the work site as the origin of coordinates, the DEM can be divided into eight direction lines and eight sectors. Given two matrices with the same row and column numbers as the DEM, they are called auxiliary grids and visual matrices, respectively. The matrix is initialized; Process the grid points on the eight direction lines and judge their occlusion relationship; The grid points within the eight sectors are processed to complete the visual field analysis of the DEM. 3. the GNSS satellite visibility forecasting method under a kind of mountainous area blocking environment as claimed in claim 2, is characterized in that, the grid point blocking relation calculation step on eight direction lines comprises: Calculate the minimum elevation value Z that is just visible to the grid point and the station; Assign values to the auxiliary grid and the visible matrix. If the grid point elevation is greater than the minimum elevation value Z, the station and grid point are visible, and the corresponding position of the grid point in the auxiliary grid and visible matrix is set to true, otherwise it is not possible. sight and false; Move the grid point outward along the direction line, repeat the above steps, if the distance reaches the buffer distance, enter the processing of the next direction line, until the calculation of the eight direction lines is completed. 4. the GNSS satellite visibility forecasting method under a kind of mountainous area shielding environment as claimed in claim 3 is characterized in that, after the direction line grid point is processed, the remaining grid points that fall in the sector are processed, and in the sector The grid point processing steps include: Calculate the minimum elevation value Z that is just visible to the grid point and the station; Assign values to the auxiliary grid and the visible matrix. If the grid point elevation is greater than the minimum elevation value Z, the station and grid point are visible, and the corresponding position of the grid point in the auxiliary grid and visible matrix is set to true, otherwise it is not possible. sight and false; Move the grid point away from the station. If the distance reaches the buffer distance, the next sector will be processed until the calculation of all eight sectors is completed, that is, the DEM visual field analysis is completed. 5. the GNSS satellite visibility forecasting method under a kind of mountainous area shielding environment as claimed in claim 1, is characterized in that, to mountain obstacle maximum elevation angle extraction and carry out polynomial fitting; Among them, the specific steps of extracting the maximum elevation angle of the mountain obstacle include: Judging any position in the visible matrix, if there are invisible points in eight adjacent positions around it, record the position. After searching all positions in the matrix, a position set is obtained, which is the boundary position of the visual field. ; Calculate the elevation angle of the mountain at the corresponding position in the position set: traverse the set, calculate the azimuth and elevation angles of the coordinates of the position in the DEM data and the coordinates of the work site jointly, and finally obtain the azimuth and elevation angles of all boundary points, denoted as { (A1,E1),(A2,E2),…,(An,En)}; Among them, the process of performing polynomial fitting includes: After extracting the maximum elevation angle of the obstacles on the mountain, multiple discrete points (A, E) can be obtained, and polynomial fitting is performed on them. The fitting criterion adopts the least squares criterion, and finally the functional relationship between the azimuth angle and the maximum elevation angle of the obstacle is obtained E = f(A), through any azimuth angle, the maximum elevation angle of the mountain obstacle at this azimuth angle can be obtained. The calculation formula includes:

In the formula, n represents the fitting order,

Represents the polynomial fit coefficients. 6. the GNSS satellite visibility forecasting method under a kind of mountainous area shielding environment as claimed in claim 1 is characterized in that, the forecast period is traversed by the set time interval deltaT, and the satellite visibility state of this moment is calculated , the processing flow of each traversal is the same. 7. the GNSS satellite visibility forecasting method under a kind of mountainous area blocking environment as claimed in claim 6, is characterized in that, described satellite visibility state calculation step is as follows: Calculate the coordinates of the satellite at time t (x ^s , y ^s , z ^s ); Calculate the altitude angle and azimuth angle of each satellite line of sight according to the coordinates of the operation site (x, y, z). The specific calculation formula is as follows:

elevation=arcsin(u) In the above formula,

Represents the unit vector from the operating site to the satellite line of sight; r represents the geometric distance from the operating site to the satellite; OMGE represents the angular velocity of the earth, which is a constant, and its value is 7.2921151467×10 ^-5 ; C is the speed of light;

Represents the vector of the station center coordinate system centered on the operating site; azimuth and elevation represent the azimuth and elevation angles of the satellites, respectively. 8. the GNSS satellite visibility forecasting method under a kind of mountainous area shielding environment as claimed in claim 7, is characterized in that, the azimuth angle of each satellite at time t is substituted into the mountain obstacle maximum elevation angle and its function with azimuth angle , the maximum elevation angle of the obstacle at time t of each satellite can be obtained; Compare the altitude angle of each satellite at time t with the calculated maximum elevation angle of the obstacle. If the altitude angle of the satellite is less than the maximum elevation angle of the obstacle, the satellite is invisible at time t; on the contrary, if the altitude angle of the satellite is greater than or equal to the maximum elevation angle of the obstacle, then This satellite is visible at time t; After the visibility of all satellites at time t is judged, the number of all visible satellites at time t is counted. 9. the GNSS satellite visibility forecasting method under a kind of mountainous area shielding environment as claimed in claim 8, is characterized in that, calculates the geometric precision factor of time t, and concrete calculation formula comprises:

Among them, n is the total number of visible satellites at time t, el in el _n represents the altitude angle, az in az _n represents the azimuth angle, n represents the nth satellite, such as el ₁ represents the altitude angle of the first satellite, az ₁ represents the azimuth of the first satellite; cos and sin represent cosine and sine functions respectively; g in the G matrix represents its element, where g ₁₁ -g ₄₄ represents its position in matrix G; GDOP, PDOP, HDOP, VDOP represents geometric precision factor, position precision factor, horizontal precision factor and vertical precision factor, respectively. 10. the GNSS satellite visibility forecast method under a kind of mountainous area blocking environment as claimed in claim 1, is characterized in that, described satellite visibility forecast result comprises the visible satellite sky view of drawing, satellite visibility time series Graph with precision factor graph.

Images