Utilizing 3-D topographical information for the quality assessment of a satellite surveying

Abstract

Fostered by the rapid development of spatial information techniques, three-dimensional topographical information (e.g., digital terrain models and digital surface models) of a higher quality is becoming accessible to the public. The main emphasis of this study is on developing a satellite visibility and quality assessment technique for a global navigation and satellite system utilizing high-resolution 3-D topographical information. An adaptive sampling and analysis procedure is proposed to increase the computational efficiency while processing high-resolution topographical data. Additionally, a unified least-squares approach is introduced to realistically model the uncertainties in satellite orbits. Finally, simulation tests and global positioning system (GPS) field work have been performed to demonstrate the performance and capability of the proposed approach. The results reveal a significant improvement in efficiency and reliability for the quality estimation of a satellite surveying. Consequently, the proposed approach will benefit applications in which a pre-evaluation of the positioning quality is of a major concern for a satellite surveying (e.g., GPS field planning or network design).

Introduction

Since the launch of its first satellite in 1978, the global positioning system (GPS) has become one of the most exciting technologies in the past decades and subsequent years. It successfully fulfills the high-quality positioning and timing needs not only in military applications as it was originally designed for but also in a wide variety of civilian uses (Leick 2004). In addition to the GPS system, more and more satellite systems which provide autonomous geospatial positioning with global or regional coverage have also been developed or are under development, for instance, the Russian GLONASS, the European Galileo, the Chinese Beidou-2/Compass, and the Indian IRNSS (see Dinwiddy et al. 2004; Leick 2004; FSARF 2005; Grelier et al. 2007; Tsai et al. 2008). The growing number of global navigation and satellite system (GNSS) satellites in space not only provides more choices to potential users but also enables better global coverage and essentially improved positioning quality.

Based on the principle of satellite positioning, the success of a GNSS surveying relies on a good network geometry constituted by satellites that are visible to a receiver. Consequently, a visibility analysis becomes an essential step in evaluating the quality of satellite positioning. The satellite visibility can be typically determined by computing the elevation angles of satellites with respect to a receiver. In a simplified scheme, the receiver is assumed to be located on an infinite plane, and satellites with positive elevation angles are regarded as visible to this receiver. A mask angle can also be imposed to exclude satellites close to the horizon. Once all visible satellites are identified, the positioning quality can be estimated using the dilution of precision (DOP) factors extracted from the covariance matrix in a least-squares adjustment (see, e.g., Strang and Borre 1997; Hofmann-Wellenhof et al. 2001; Misra and Enge 2001; Leick 2004). This approach usually gives a satisfactory quality of evaluation in open areas. However, in a case which involves complicated surface variations (e.g., an urban or mountain area), the quality is usually overestimated due to the fact that many satellite signals are actually obstructed by buildings and/or mountains (Li and Han 2008). The inclusion of realistic surface information in the satellite positioning analysis is thus essential, especially for applications in an urban or mountain area. Recently, more and more studies have focused their attention on topographical and obstruction issues. For example, Taylor et al. (2005) utilized a triangular irregular network created by LiDAR and photogrammetry data to model the buildings while analyzing possible signal obstructions. Hogan and Santos (2005) studied the data link between the GNSS base and rover stations with a digital terrain model. Li et al. (2006) implemented a ray tracing model to investigate the GPS multipath effect using a high-resolution digital surface model (DSM). Furthermore, Zhang et al. (2008) and Kleijer et al. (2008) used high-fidelity city models to investigate GNSS availability and accuracy in urban areas. These preceding works, although having various purposes, have demonstrated the advantages of incorporating digital topographical information in a satellite survey. However, the inclusion of detailed surface information may also degrade the computational efficiency which becomes another issue needed to be further discussed.

In addition to the satellite visibility and network configuration, the uncertainty in satellite orbits is also a key factor that affects positioning quality. Taking the GPS for an example, its satellite orbit solutions provided by the International GNSS Service (IGS) come with five different levels of quality. The uncertainties range from 2.5 to 100 cm in satellite positions and from 75 ps to 5 ns in satellite clocks (IGS 2010, see Table 1). If orbits with different levels of uncertainty are used for estimating the receiver’s positions, one can expect solutions to be of variable quality. Consequently, when the quality of GNSS positioning is being assessed, the orbit uncertainties should also be properly considered in order to yield a realistic result. This can be done by treating the satellite orbits as observables so that their uncertainties are included during the adjustment computations. A unified approach to the least-squares adjustment can be applied to properly reflect the uncertainties from all variables. This approach treats all variables in the mathematical model as observations with random errors such that all associated uncertainties can be realistically included in the adjustment (Mikhail and Ackermann 1976). Owing to the flexibility of this approach, the satellite orbits can provide a partial datum constraint (i.e., soft constraint) at a level based on the orbit uncertainties, and thus a more realistic analysis can be performed.

Table 1 Different GPS orbits provided by IGS (2010)

Satellite visibility analysis

In order to perform a satellite visibility analysis, the positions of satellites should be first estimated using almanac or ephemeris files. The detailed algorithm for computing satellite positions based on the normal orbit theory can be found in Leick (2004) which is not repeated here. With known satellite positions, the visibility of satellites can be determined by a line-of-sight (LOS) analysis between a receiver and the satellites using three-dimensional topographical data.

The basic method of a LOS analysis is to determine whether any obstruction blocks the sight vector between a view point (i.e., receiver) and a target (Guth 2004), as illustrated in Fig. 1. The first step of this method is to slice the local terrain into pieces in every direction from the receiver (Fig. 2). Next, points are sampled on each slice, and all the elevation angles of each sampling point with respect to a receiver are calculated. Then, the obstruction angle in each direction is determined by finding the maximum elevation angle on each slice. As shown in Fig. 3, a vector between a receiver and a target point with an elevation less than the obstruction angle will be blocked by the terrain. Consequently, the target is not visible to that receiver.

Fig. 1

LOS analysis on a terrain profile

Fig. 2

Sampling slices around the receiver for LOS analysis

Fig. 3

The maximum elevation angle on each slice for determining satellite visibilities

By repeating the above process for all directions around the receiver, the maximum elevation angle in each azimuth is obtained and can be plotted as a closed curve line in a sky plot. Figure 4 shows the obstruction line in a sample sky plot. Every satellite position can be checked to see if it is located inside (visible) or outside (invisible) the obstruction line, and the satellite visibility can be determined accordingly.

Fig. 4

The obstruction (black line) in a sky plot

In practice, this visibility analysis algorithm is time-consuming, especially when digital terrain data of large sizes are involved. In order to improve the efficiency of the analysis, an adaptive sampling algorithm for the satellite visibility analysis is proposed and described in the next section.

Adaptive topographical analysis algorithm

The sampled points on the terrain provide basic information for a visibility analysis. It is easy to realize that a higher sampling resolution should give a better result in the analysis. However, it could become inefficient when analyzing digital terrain data sets of large sizes using a high sampling resolution. An improved sampling algorithm should be developed to increase computational efficiency while maintaining acceptable quality in the result of the analysis.

One of the important parameters to determine the sampling interval is the orbit resolution angle from a receiver. Its value is proportional to the sampling interval of a satellite orbit and inversely proportional to the distance between the receiver and the satellite (see Fig. 5). For a typical GPS orbit with a 5-min sampling interval, the minimum orbit resolution angle (θ _s ) is about 3.3° for any receiver on the earth’s surface. Furthermore, when ground points are analyzed with a certain sampling interval, it produces a sampling resolution angle (θ) from a receiver to sampled points. It is illustrated in Fig. 6 that this angle is dependent on the sampling interval d _s , the distance d, the elevation angle El _A , and the slope S _hl,h2 at the sampled point. The mathematical expression has been derived to relate these variables and can be written as Eq. 3.1. Apparently, a proper choice of the sampling interval should result in a sampling resolution angle which is equal or less than the satellite orbit resolution angle (i.e., θ ≤ θ _s ).

Fig. 5

The viewing angle of a receiver with respect to two satellite sampling epochs

Fig. 6

a Graphical relationships between all parameters required for determining the sampling interval and b a closed view showing detailed geometry of these parameters

$$ {d_s} = \frac{{d \times \sin \theta \times \cos {S_{h1.h2}}}}{{\sin \left( {{S_{h1.h2}} - E{l_A} - \theta } \right) \times \cos E{l_A}}} $$

(3.1)

In a practical application with a given satellite orbit, a proper viewing angle θ can be estimated and is thus a known constant. According to Eq. 3.1, the sampling interval does not stay constant but depends on the distance of the sampling point to the receiver (d) and the terrain complexity at each sampling point (characterized by El _A and S _hl,h2), as shown in Fig. 7.

Fig. 7

Varied intervals resulting from the adaptive sampling algorithm

Positioning quality assessment with orbital uncertainties

Two important factors have been considered in evaluating the quality of a satellite surveying. They are the number of visible satellites and the geometry constituted by these visible satellites with respect to a receiver. Nevertheless, a satellite orbit always comes with a certain level of uncertainty. It should be included in the quality assessment in order to give a more realistic result.

Figure 8 depicts a basic form of GPS single-point positioning. The satellites $({X^{si}},{Y^{si}},{Z^{si}})$ are assumed to be known points in the sky, and the receiver $({X_r},{Y_r},{Z_r})$ is an unknown point on the ground. The location of the receiver $({X_r},{Y_r},{Z_r})$ can be solved by the intersection of range observations P _i (i = 1∼ n) between the satellites and the receiver. When the clock offset dt is considered, one can write an observation equation for each range measurement as:

Fig. 8

Single-point GPS positioning with range measurements

$$ {P_i} - \sqrt {{{{\left( {{X^{si}} - {X_r}} \right)}^2} + {{\left( {{Y^{si}} - {Y_r}} \right)}^2} + {{\left( {{Z^{si}} - {Z_r}} \right)}^2}}} + c \cdot dt = 0\quad c:{\hbox{speed}}\,{\hbox{of}}\,{\hbox{signal}} $$

(4.1)

According to Eq. 4.1, the satellite position can be uniquely determined if four range measurements are available. In the case when n satellites are visible to a receiver, one can write the following equations:

$$ \left\{ {\begin{array}{*{20}{c}} {{F_1} = {P_1} - \sqrt {{{{\left( {{X^{s1}} - {X_r}} \right)}^2} + {{\left( {{Y^{s1}} - {Y_r}} \right)}^2} + {{\left( {{Z^{s1}} - {Z_r}} \right)}^2}}} + c \cdot dt = 0} \\{{F_2} = {P_2} - \sqrt {{{{\left( {{X^{s2}} - {X_r}} \right)}^2} + {{\left( {{Y^{s2}} - {Y_r}} \right)}^2} + {{\left( {{Z^{s2}} - {Z_r}} \right)}^2}}} + c \cdot dt = 0} \\{{F_3} = {P_3} - \sqrt {{{{\left( {{X^{s3}} - {X_r}} \right)}^2} + {{\left( {{Y^{s3}} - {Y_r}} \right)}^2} + {{\left( {{Z^{s3}} - {Z_r}} \right)}^2}}} + c \cdot dt = 0} \\{{F_n} = {P_n} - \sqrt {{{{\left( {{X^{sn}} - {X_r}} \right)}^2} + {{\left( {{Y^{sn}} - {Y_r}} \right)}^2} + {{\left( {{Z^{sn}} - {Z_r}} \right)}^2}}} + c \cdot dt = 0} \\\end{array} } \right. $$

(4.2)

which can be written symbolically as:

$$ {\mathbf{F}}\left( {{\mathbf{x}},{\mathbf{l}}} \right) = 0 $$

(4.3)

where 1 is an (n × 1) observation vector formed by all range measurements P _i , and x is an (4 × 1) unknown parameter vector containing the three coordinates of a receiver and the clock offset. Generally, this equation can be linearized as:

$$ {\mathbf{Bv + A\Delta }} = {\mathbf{f}} $$

(4.4)

where ${\mathbf{B}} = {\left. {\frac{{\partial {\mathbf{F}}}}{{\partial {\mathbf{l}}}}} \right|_{{{\mathbf{x}}_0},{{\mathbf{l}}_b}}}$and ${\mathbf{A}} = {\left. {\frac{{\partial {\mathbf{F}}}}{{\partial {\mathbf{x}}}}} \right|_{{{\mathbf{x}}_0},{{\mathbf{l}}_b}}}$ are coefficient matrices with respect to the observations and unknown parameters; v is a residual vector for the observations 1 _b ; ∆ is a correction vector for approximate parameter values x ₀, and f is a constant vector. If all usable satellite positions are treated as observables (denoted by ${{\mathbf{l}}_{{{\mathbf{x}}_S}}}$ with a residual${{\mathbf{v}}_{{{\mathbf{x}}_S}}}$), the observation vector is extended as:

$$ {\mathbf{\bar{l} = }}{\left[ {\begin{array}{*{20}{c}} {{{\mathbf{l}}^{\mathbf{t}}}} & {{\mathbf{l}}_{{{\mathbf{x}}_{\mathbf{S}}}}^{\mathbf{t}}} \\\end{array} } \right]^t} $$

(4.5)

Consequently, Eq. 4.4 should be rewritten as:

$$ {\mathbf{\bar{B}\bar{v} + A\Delta }} = {\mathbf{f}} $$

(4.6)

where ${\mathbf{\bar{B}}} = {\left. {\frac{{\partial {\mathbf{F}}}}{{\partial {\mathbf{\bar{l}}}}}} \right|_{{{\mathbf{x}}_0},{{{\mathbf{\bar{l}}}}_b}}}$ and ${\mathbf{\bar{v} = }}{\left[ {\begin{array}{*{20}{c}} {{{\mathbf{v}}^{\mathbf{t}}}} & {{\mathbf{v}}_{{{\mathbf{x}}_{\mathbf{S}}}}^{\mathbf{t}}} \\\end{array} } \right]^{\mathbf{t}}}$. The corresponding a priori cofactor matrix for the observables is:

$$ {\mathbf{\bar{Q} = }}\left[ {\begin{array}{*{20}{c}} {\mathbf{Q}} & 0 \\0 & {{{\mathbf{Q}}_{{{\mathbf{x}}_{\mathbf{S}}}{{\mathbf{x}}_{\mathbf{S}}}}}} \\\end{array} } \right] $$

(4.7)

where Q is the a priori cofactor matrix for range measurements and ${{\mathbf{Q}}_{{{\mathbf{x}}_{\mathbf{S}}}{{\mathbf{x}}_{\mathbf{S}}}}}$ is the a priori cofactor matrix for all usable satellite positions. The unified least-squares solution to this model can thus be obtained by (Mikhail and Ackermann 1976):

$$ {\mathbf{\Delta }} = {\left( {{\mathbf{N + Q}}_{{\mathbf{xx}}}^{ - {\mathbf{1}}}} \right)^{ - 1}}{\mathbf{\dot{t}}} $$

(4.8)

where ${\mathbf{N = }}{{\mathbf{A}}^{\mathbf{t}}}{\mathbf{Q}}_e^{ - 1}{\mathbf{A}}$, ${\mathbf{t = }}{{\mathbf{A}}^{\mathbf{t}}}{\mathbf{Q}}_e^{ - 1}{\mathbf{f}}$, ${{\mathbf{Q}}_{\mathbf{e}}} = {\mathbf{\bar{B}\bar{Q}}}{{\mathbf{\bar{B}}}^{\mathbf{t}}}$, ${\mathbf{\dot{t} = }}{{\mathbf{A}}^{\mathbf{t}}}{\mathbf{Q}}_{\mathbf{e}}^{ - {\mathbf{1}}}{\mathbf{f + Q}}_{{\mathbf{xx}}}^{ - {\mathbf{1}}}{{\mathbf{x}}_0}$, and Q _xx is the a priori cofactor matrix of the parameters to be solved. Finally, the adjusted parameters can be computed by:

$$ {\mathbf{\hat{x} = }}{{\mathbf{x}}_{\mathbf{0}}}{\mathbf{+ \Delta }} $$

(4.9)

with their cofactor matrix written as:

$$ {{\mathbf{Q}}_{{\mathbf{\hat{x}\hat{x}}}}} = {\left( {{\mathbf{N + Q}}_{{\mathbf{xx}}}^{ - {\mathbf{1}}}} \right)^{ - 1}} $$

(4.10)

As illustrated in the above formulations, the quality of the final solution depends on the a priori uncertainties of three sets of variables: the range measurements (denoted by Q), the satellite orbits (denoted by ${{\mathbf{Q}}_{{{\mathbf{x}}_{\mathbf{S}}}{{\mathbf{x}}_{\mathbf{S}}}}}$), and the unknown parameters (denoted by Q _xx ). In other words, the complete set of uncertainties from all variables in the model has been realistically taken into account which will help to achieve a more reliable estimate of accuracy for the adjusted parameters.

This cofactor matrix is usually transformed to a local geodetic coordinate system so that the DOP factors can be estimated for a receiver’s horizontal and vertical positions and its corresponding time precision.

Numerical validations

Performance test of the adaptive sampling algorithm

In this test, the proposed adaptive sampling algorithm is applied to process a real DSM data set covering 3,800 by 4,100 m with a 2-m resolution (see Fig. 9). The sampled points are used for a visibility analysis, and the maximum viewing angle resolution is set to be 1°. The same test has also been performed for a regular sampling algorithm (sampling interval = 2 m). The results from both approaches are represented as obstruction lines in sky plots (see Fig. 10). A quantitative comparison is also listed in Table 2.

Fig. 9

The DSM data used for a performance test of the proposed adaptive algorithm

Fig. 10

The sky plots with obstruction lines created from a a regular sampling algorithm and b the proposed adaptive sampling algorithm

Table 2 Performance comparison of the regular and adaptive sampling algorithms

From Fig. 10, the predicted obstructions in a sky plot from the regular and adaptive sampling algorithms are visually identical. Furthermore, it is shown in Table 2 that the required processing time for a visibility analysis reduces significantly from 99.735 to 10.870 s when the adaptive algorithm is applied. It is also illustrated that the root mean square difference (RMSD) between the estimated obstruction elevation angles from the two algorithms is 0.05°, which is insignificant compared to the preset 1° viewing angle resolution. In other words, the proposed algorithm has demonstrated its capability of improving computational efficiency while keeping the visibility analysis to an acceptable quality.

Reliability test for the positioning quality assessment model

In order to check whether the proposed quality assessment approach gives reasonable solutions, a simulation test was performed. The satellite orbits were created from almanac data (YUMA 449). The range measurements from satellites to a receiver located at (λ = 121.56443⁰, φ = 25.03367⁰) were also simulated. Two-centimeter random errors were imposed on the range measurements, and 0.01- to 0.1-m random errors were added to the satellite orbits before they were used to estimate the receiver’s coordinates. The estimated receiver coordinates were compared to their true (simulated) values, and the total root mean square errors were plotted as a red line in Fig. 11. In the same figure, the blue line represents the predicted accuracies of the receiver’s coordinates by the proposed approach, and the green line denotes the predicted accuracies of the receiver’s coordinates by an approach without taking into account the orbital uncertainties.

Fig. 11

True errors vs. predicted errors in single-point positioning

The results show that the accuracies predicted by the proposed approach exhibit a reasonable magnitude and trend compared to the true errors. On the other hand, the comparison approach does not give a realistic accuracy estimate since the orbital uncertainties were not included in the analysis.

A field test

The test area is located in the Wenshan area of Taipei city in Taiwan (bottom left at λ = 121.57249°, φ =24.97292°; top right at λ =121.60331°, φ =24.99808°). Two sets of DSM data covering this area are available. One is a 5-m resolution DSM in which most buildings are artificially removed, and the other DSM is generated from LiDAR scan data with a 2-m resolution. Figure 12 displays these two sets of DSM in a grayscale image.

Fig. 12

a 5 m DSM and b 2 m LiDAR DSM for the test area shown in grayscale maps

Satellite visibilities for a 24-h period on April 1, 2009, have been analyzed for six test sites in this area using the proposed approach and the GPS almanac file at week 502. Figure 13 shows the site locations and predicted numbers of visible satellites for these six sites. The cyan lines and blues lines show the results using the 5-m DSM and 2-m DSM, respectively. The red lines present the results without any topographical consideration (assuming a flat plane).

Fig. 13

Locations and numbers of visible satellites for the six test sites

As shown in Fig. 13, the visibility analysis results using DSM data are significantly different from those obtained without a topographical consideration. In each test site, the number of visible satellites decreases when DSM data are included in the analysis, indicating that topographical variations have a nonnegligible impact on the satellite visibilities. Additionally, the number of visible satellites using the high-resolution (2-m) DSM is less than that using low-resolution DSM. It is due to the fact that most buildings have been removed in the low-resolution (5-m) DSM. Consequently, the obstruction effect is underestimated when the low-resolution DSM is used in the analysis.

To further verify the above analysis results, a GPS field survey was carried out simultaneously at site 5 (with a higher topographical obstruction) and site 6 (with a lower topographical obstruction) for a 2.5-h period from 10:30 a.m. to 12:00 p.m. on April 1, 2009, local time. The observation data were converted to the receiver independent exchange format, and the visibility was retrieved and plotted together with the above-mentioned results in Figs. 14 and 15.

Fig. 14

a Visibility test result for site 5 and b its close view during the period of the field survey

Fig. 15

a Visibility test result for site 6 and b its close view during the period of the field survey

From both Figs. 14 and 15, it is clear that the actual satellite visibilities at these two sites show good agreement with the predictions from DSM data but were significantly different to those without topographical considerations. The results were further analyzed quantitatively using a rate of prediction index which is defined as:

$$ {I_{PD}} = \frac{{{n_o}}}{{{n_p}}} $$

(5.1)

where n _o denotes the actual number of observable satellites and n _p represents the number of observable satellites predicted from a visibility analysis. The comparison results are listed in Table 3.

Table 3 The rate of prediction in the field test

In Table 3, the satellite visibilities are enormously overestimated when the topography (DSM) is not included in the analysis. At site 6, where the actual obstruction is lower, the rate of overestimation is 27%. At site 5 which is located in a higher-obstruction area, the rate of overestimation reaches 103%. These overestimations are significantly reduced when the DSM is included in the analysis. Using the low-resolution DSM reduces the overestimation to 30% for site 5 and to 9% for site 6. Using the high-resolution DSM further reduces the overestimation to be negative 8% and 7% at these two sites (underestimation). In summary, using the high-resolution DSM gives the best visibility analysis result for both high-obstruction and low-obstruction areas. The low-resolution DSM also gives a good result in a lower-obstruction area but becomes less realistic in a higher-obstruction area.

Finally, the sky plots of the predicted and observed visibilities were generated and shown in Fig. 16. These sky plots again verify that the proposed approach produces a realistic result for satellite visibilities.

Fig. 16

a The predicted sky plots and b the observed sky plots of site 5 (left) and site 6 (right)

Conclusions

Recent developments in the GNSS technique make it more popular and significant in people’s daily life. A fast and accurate positioning solution is actively pursued in both engineering applications and scientific research. In this study, we have demonstrated the feasibility of using high-fidelity digital topographical information in the GNSS satellite positioning. According to the proposed approach, the positioning quality can be realistically and efficiently assessed. Consequently, a field schedule of a satellite surveying project can be precisely planned in the office, ensuring that a GNSS surveying of good quality is to be performed. Furthermore, although this study focuses on a single-point positioning scenario, the proposed approach can be directly applied to the case of relative positioning. Relative positioning requires simultaneous visibilities from multiple receivers to common satellites, such that the obstruction effect due to terrain variations and its impact on positioning quality is much more complicated. With the help of high-fidelity digital topographical data and the proposed analysis approach, the quality can always be visualized clearly and precisely for a GNSS surveying.