Unified Astronomical Positioning and Orientation Model Based on Robust Estimation

Astronomical geodesy is a technique that uses the sun, moon, stars, and other natural objects as beacons; takes the celestial horizon coordinates of such objects as observations; and finally determines the geographical position of the station and the azimuth of one direction (Hirt et al. 2010). Because it is a kind of absolute positioning, astrogeodetic points in the geodetic control network are often used to control the accumulation of measurement error and determine vertical deflection (Hirt and Flury 2008; Hirt and Seeber 2008; Müller et al. 2004). Because of their advantages, such as good stealth, not easy to be disturbed, and strong independence, astronomical positioning and orientation could be used as a method of independent autonomous navigation (Kaplan 1999; Vulfovich and Fogilev 2010). In the traditional method, astronomical positioning and orientation are realized by theodolite to observe natural objects one by one; the horizontal angle and elevation angle are recorded at the same time; and then the astronomical longitude, latitude, and azimuth are solved, respectively, according to algorithms (Robbins 1977). In the celestial sphere, the north celestial pole, zenith, and stars can constitute a positioning triangle, as shown in Fig. 1.

In Fig. 1, $P$ = north celestial pole; $Z$ = zenith; $W$ = west point; $WMQ$ = horizon circle; $WNQ\prime$ = celestial equator; and $\sigma$ = a star in the celestial sphere. According to the spherical triangle cosine formula, the relationship among zenith angle $z$, astronomical latitude $\phi$, longitude $\lambda$, star right ascension $\alpha$, declination $\delta$, and Greenwich true sidereal time $S$ iswhere, the zenith distance $z$ need to be corrected by atmospheric refraction according to a well-known equation (Meeus 1998). The Greenwich true sidereal time $S$ could be converted from observing epoch $T_{\text{UTC}}$. Traditionally, $T_{\text{UTC}}$ is obtained by radio timing, which needs to receive the radio time signal broadcasted by the timing center. Currently, it can be obtained by synchronizing the measuring pulse to a global navigation satellite system (GNSS) receiver, and then the observing epoch is fixed on the Universal Time Coordinated (UTC) scale.

In Eq. (1), right ascension $\alpha$ and declination $\delta$ can be obtained through an apparent position calculation procedure (Bangert et al. 2009). Zenith distance $z$ and UTC time $T_{\text{UTC}}$ can be obtained through observation. Therefore, only $\lambda$ and $\phi$ remain unknown in Eq. (1); if $n$ stars have been observed, the astronomical longitude and latitude can be calculated according to the least-squares method, and then astronomical positioning can be realized.

When the astronomical longitude and latitude of the observing location are known, there are several kinds of astronomical orientation methods to determine the azimuth of a ground target, but the most frequently used one is the Polaris hour angle method (Adams 1968).

As is shown in Fig. 2, the azimuth ($A_B$) of ground target $B$ is the angle between the meridian plane and the vertical plane that goes through the station and the target two points; it is measured clockwise from north point $N$ in the range of 0–360°. Although north point $N$ cannot be observed directly like ground targets, it can be determined by means of observing the azimuth of celestial bodies (Zhan et al. 2016; Carter et al. 1978). Let $L_{\sigma }$ represent the observed horizontal angle of celestial body at Greenwich UTC time $T_{\text{UTC}}$, let $L_B$ represent the observed horizontal angle of the ground object, let $A_{\sigma }$ represent the azimuth of Polaris at the observing epoch, and let $L_N$ represent the horizontal angle of north point $N$; therefore:

Obviously, the azimuth $A_B$ of the ground target can be determined based on the azimuth $A_{\sigma }$ of Polaris at the observing epoch. According to the principle of hour angle method, $A_{\sigma }$ can be calculated according to the following formula:

Then the azimuth $A_B$ can be calculated according to Eq. (2). In practice, the final azimuth can be estimated as the average of $A_B$ if the Polaris has been observed $m$ times, and then astronomical orientation is realized.

With the development of imaging and digital processing technology, a new kind of celestial navigation method that uses star sensors on the ground and at sea to image multiple stars simultaneously has appeared (Hughes et al. 2009; Li et al. 2014; Parish et al. 2010). Unfortunately, some defects still exist in these methods. On the one hand, astronomical positioning and orientation are implemented separately, leading to low efficiency. On the other hand, although multiple stars can be imaged simultaneously by a charge-coupled device (CCD), some low-precision star observations still exist because of weather or other severe conditions, and calculating with these data at the same weight directly will reduce the accuracy of astronomical positioning and orientation.

To solve these problems, a unified astronomical positioning and orientation model based on robust estimation is proposed. By using the Rodrigues matrix, this model realizes the calculation of astronomical longitude, latitude, and azimuth simultaneously. By introducing robust estimation, the influence of outliers is excluded, the use of low-precision observations is limited, and the high-precision observations are fully utilized; thus, more reliable and effective valuations of astronomical longitude, latitude, and azimuth can be obtained.

In the traditional method, astronomical positioning is generally carried out as the first step, and then the calculated longitude and latitude are deemed as known quantities to determine the astronomical azimuth of a ground target; thus, errors in longitude and latitude will further increase the error of the azimuth. In the horizontal circle astronomy method, the error equation is established on the difference between the observed and calculated azimuth, and the longitude, latitude, and azimuth can be solved simultaneously (White 1966). However, the elevation-angle observations are underutilized. Therefore, the horizontal-angle and elevation-angle observations are both used in the unified astronomical positioning and orientation model. This method is actually a vector algorithm that can not only improve the calculation efficiency and uniformity but also restrict the error propagation and improve the calculation accuracy.

Suppose at the observing epoch $T_{\text{UTC}}$, the horizontal angle and zenith distance of a celestial body are, respectively, $L$ and $z$. The rectangular coordinate of the celestial body in the instrument coordinate system is

In the true equatorial coordinate system, the coordinates of the celestial body could be represented by right ascension $\alpha$ and declination $\delta$ aswhere $(\alpha ,\delta )$ is the apparent position of the celestial body, which can be calculated according to the star catalog, and it can be regarded as a known quantity. Therefore, the rotational transform relationship between $C_1$ and $C_2$ can be expressed as (Capon 1954)where $A_H$ = azimuth of the $X_H$ axis in the horizontal coordinate, namely, the azimuth of horizontal-angle value 0; $P_y$ = reverse matrix of the $Y$ axis; and $R_y({90}^{\circ }-\phi )$ and $R_z(S+\lambda )$ represent the corresponding rotation matrix.

In Eq. (9), $S$ represents the Greenwich true sidereal time of the observing epoch, which can be converted from the UTC time $T_{\text{UTC}}$. Therefore, Eq. (6) could be changed to the following:

If order:

Then:where the specific form of $R$ can be expanded as

Obviously, it is difficult to calculate $\lambda$, $\phi$, and $A_H$ according to Eqs. (13) and (14) directly. The antisymmetric matrix $Q$ is introduced here as follows:

Because $R$ is an orthogonal rotation matrix with three degrees of freedom, it can be seen as a Rodrigues matrix that is composed of $Q$ and $I$ (Yao et al. 2006), namely:

Then Eq. (13) is changed as follows:

It can be further expressed in the form of rectangular coordinates as

Antisymmetric matrix $Q$ has the following property:

Multiply ${(I-Q)}^{-1}$ on the right side:

Then multiply ${(I-Q)}^{-1}$ on the left side:

Substitute Eq. (21) into Eq. (18):

Then substitute Eq. (15) into Eq. (22) and expand it:

Because there are only two independent equations in Eq. (23), the three independent parameters $a$, $b$, and $c$ cannot be solved by observing only one celestial body. However, if two or more celestial bodies are observed, then four or more than four independent equations can be determined; the three independent parameters $a$, $b$, and $c$ can be calculated by the least-squares method; and then the matrix $Q$ and the rotation matrix $R$ can be obtained accordingly. Therefore, the astronomical longitude $\lambda$, astronomical latitude $\phi$, and azimuth $A_H$ can be expressed according to Eq. (14):

Then, in the horizontal plane, the azimuth $A_B$ of the observed ground target is

Thus, astronomical positioning and astronomical orientation are realized simultaneously in this unified model.

Because of the influence of severe conditions, such as clouds and background light, some errors and even few outliers inevitably exist in the observations, which will cause the final least-squares estimations of longitude, latitude, and azimuth to deviate from their true values. Considering this problem, robust estimation, which helps to take full advantage of high-precision observations, limit low-precision ones, and exclude outliers, is introduced. This way, more reliable, effective, and meaningful parameter estimations can be obtained (Li et al. 2013).

Robust estimation stems from the concept of statistical robustness, and it is proposed for least-squares estimation that has no characteristic of anti-interference. In Eq. (23), if order:then Eq. (23) can be simplified as an error equation:where $V$ = residual vector. According to the principle of robust estimation, the solution of Eq. (27) is:

Where $\overline{P}$ is the equivalent-weight matrix, the success of the robust estimation mainly depends on the accuracy of the initial equivalent weights. The observations in astronomical surveying are independent from each other; therefore, $\overline{P}$ is a diagonal matrix, and its diagonal element is ${\overline{p}}_i$ (Yang et al. 2002). Because the accuracy of each observation is not equal, and observations may even contain gross errors, the least-squares estimations are not suitable as the initial value. According to the characteristics of astronomical measurements, first set the $L_1$ norm minimum as a criterion to calculate the initial parameters (Yang et al. 2001). Suppose:where $p_i$ = a priori weight of each observation; and $v_i$ = corresponding least-squares residual. The initial parameters without the impact of outliers can be calculated by Eq. (29); the initial residual of each observation can be calculated according to Eq. (27); and the initial estimated values of astronomical longitude, latitude, and azimuth included in the vector $\widehat{X}$ can be solved by Eq. (28). However, to get the final solution, the iterative method is generally used; the $n$th iterative solution iswhere ${\overline{P}}^n$ = $n$th equivalent-weight matrix, which should be designed elaborately. Astronomical measurement practice indicates that the main part of the observations still obeys a normal distribution, and there are just a few suspicious observations or rare abnormal observations. Therefore, in the process of robust estimation, all the reliable information should be fully used and keep their original weights. Nevertheless, the weight of suspicious observations should be reduced according to their degree of credibility, and the outliers should be eliminated. From the principle of robust estimation, a famous equivalent-weight function is the IGG3 scheme (Yang 1999); the equivalent-weight equation is

Normally, in Eq. (31), the value of $k_0$ and $k_1$ can be assigned as 1.5 and 3.0, respectively; ${\tilde{v}}_i{}^{n-1}$ is the (n – 1)th standardized residual of the $i$th observation, namely:where $v_i{}^{n-1}$ = observation residual; and $m_{v_i{}^{n-1}}$ = corresponding mean square error of $v_i{}^{n-1}$, which can be calculated as

In Eq. (33), ${\sigma }_0$ is the unit-weight mean square error; its theoretical or empirical value is often used. ${\overline{P}}^{n-1}$ = (n – 1)th equivalent-weight matrix; $p_i{}^{n-1}$ = $i$th diagonal element in ${\overline{P}}^{n-1}$; and $A_i$ = $i$th line of matrix $A$—namely, it is the $i$th coefficient vector of the error equation (Yang and Wu 2006). Therefore, the final robust valuation of the parameters can be calculated by Eq. (30). The astronomical longitude $\lambda$, latitude $\phi$, and azimuth $A_B$ can be calculated according to Eqs. (24) and (25).

To verify the effectiveness and accuracy of the proposed unified astronomical positioning and orientation model based on robust estimation, two kinds of observations were used to calculate the longitude, latitude, and azimuth, respectively. The first kind of observation was obtained by a large field of view (FOV) star sensor placed on the ground, which simply needed to capture a star image and eliminate artificial operation when carrying out the experiment. After being processed by such automated steps as star extraction, determination of star image center coordinates, and star pattern matching, a total of 541 stars was identified from the star image [Fig. 3(a)]. The other kind of observation was obtained by a Leica total station (TS 50). During the experiment, 24 stars were manually observed in a period, and each star was measured approximately 10 times [Fig. 3(b)]. The two experiments were carried out at the same position; the true value of the instrument station ($\lambda {\hspace{0.17em}=113}^{\mathrm{°}}6\prime 13.50″$, $\phi =34°31\prime 28.39″$) and azimuth ($A_H=43°19\prime 56.27″$) of a well-designed ground target were determined by first-class astronomical surveying, and their accuracies were 0.3″ and 0.5″, respectively (General Administration of Quality Supervision, Inspection and Quarantine of the People’s Republic of China 2000). One of the most important details is how to determine the observing time because errors in the observing epoch will directly propagate into astronomical longitude. First, the authors used the pulses per second (1 pps) single and the National Marine Electronics Association-0183 (NMEA-0183) data generated by a timed GNSS receiver to correct a time counter in the measure control unit before starting measurement, and then the time counter was synchronized to UTC scale; therefore, the GNSS receiver is not needed in the process of measuring. In the star sensor experiment, the precise time delay of the shutter is necessary. A calibration procedure was developed based on the comparison of two epochs. The first one is a count value recorded as the shutter opening epoch by the control unit, and the other one is imaged by the star sensor from a screen displaying the real-time value of the time counter. The comparison of both count values revealed a significant delay of the shutter motion, and then the time delay could be calculated The accuracy of the implemented timing technique has been verified to be approximately ±5 ms. In the total station experiment, the observing epoch was determined by an electronic pulse from the control unit that was synchronously sent to the corrected time counter and to the total station; the accuracy of this method has been demonstrated to be approximately ±3 ms.

As can be seen from Fig. 3(a), the stars are mainly uniformly distributed throughout the sky. However, in Fig. 3(b), the stellar azimuths are different from each other, resulting in the elevation angle of each star appearing as a column. It can be seen from Fig. 3(b) that the elevation angle of each star is approximately 60°, and the horizontal-angle intervals are similar to each other, indicating that the 24 stars are substantially distributed near the 60° altitude circle.

According to the star number, the theoretical azimuth and elevation angles of these stars at their observing epoch can be calculated accurately. The observation azimuth of the stars can be obtained by adding the observation horizontal angles on the basis of the azimuth of the star sensor’s x-axis, and then by subtracting it from the true azimuth of these stars, the residual of the horizontal-angle observations can be obtained; these are shown in Fig. 4.

As can be seen in Fig. 4, although the observed precision of the total station was significantly higher than that of the large FOV star sensor, they both obeyed normal distribution. The horizontal angle residuals of most stars in Fig. 4(a) were found to be within the range of ±100″, the observed accuracy of which is relatively high, and a few of the azimuth errors were within the range of ±100–±300″, the observed accuracy of which is relatively low. Similarly, the horizontal-angle residuals in Fig. 4(b) were found to be mostly within the range of ±4″, and the others were within the range of ±4–±10″.

The elevation angles of the two kinds of observations were both influenced by atmospheric refraction, so they should be corrected before being used to calculate the parameters. However, the elevation angles cannot be completely corrected by the refraction correction equation; it always remains a quantity in the elevation-angle values. Therefore, the residuals of elevation angles from their true values may still contain some errors, but their regularities of distribution can be used to research the data quality (Fig. 5).

As shown in Figs. 5(a and b), the elevation-angle residuals of most stars were found to be within the range of ±50″ to ±5″, respectively; a few of the elevation angle residuals were within the range of ±50″ to ±110″ and ±5″ to ±10″, respectively; and individual errors were greater than ±110″ and ± 10″, respectively, which could be considered as outliers.

Three processing schemes were used to deal with the experimental data. Scheme 1 used the least-squares (LS) algorithm to calculate the astronomical longitude and latitude first, then used the hour angle method to determine the azimuth. Scheme 2 used the unified method based on the Rodrigues matrix to solve the astronomical longitude, latitude, and azimuth simultaneously. Scheme 3 further adopted robust estimation based on Scheme 2 to calculate the parameters. The difference between the results of the three schemes and their true values are shown in Table 1.

As can be seen from Table 1, the positioning error and orientation error of Scheme 2 were close to those of Scheme 1, which indicates that the proposed model using the Rodrigues matrix in the solving process can realize astronomical positioning and orientation simultaneously. Furthermore, by comparing the positioning error and orientation error of Scheme 3 with those of Scheme 2, it can be determined that the positioning error declined to 2.48″ from 4.27″, and the orientation error declined to –5.64″ from –8.13″ for the star sensor observations. This is because 17 observations were identified as outliers and were given zero weight; 164 other observations were identified as abnormal values, and their weights were reduced; and another 360 high-precision observations were given greater weight in the robust estimation process. Similarly, 9 total station observations were identified as outliers, 49 other observations were identified as abnormal values, and another 204 observations were identified as high-precision values. Therefore, the astronomical positioning error declined to 0.50″ from 1.06″, and the orientation error changed to 2.05″ from 3.37″ for the total station observations. These results show that the proposed unified model based on robust estimation can effectively take full advantage of high-precision observations, limit utilization of low-precision observations, and exclude outliers, through which more accurate astronomical longitude, latitude, and azimuth estimations can ultimately be obtained.

A unified astronomical positioning and orientation model based on robust estimation is proposed. Through the introduction of the Rodrigues matrix, this model achieves the calculation of astronomical longitude, latitude, and azimuth simultaneously, eliminating the inconvenience of the traditional method in which astronomical positioning and orientation must be carried out separately in. In addition, by the use of a robust estimation method, this model can effectively suppress the influence of outliers, reduce the impact of systematic errors, make full use of high-precision observations, and eventually obtain results with relatively high accuracy. Compared with traditional models, the proposed method has improved efficiency and accuracy. Even in foggy or cloudy weather conditions, the astronomical positioning and orientation problems can be solved simultaneously by just observing several visible bright stars; therefore, the method also has important application prospects in the celestial navigation field.

This research is funded by the National Nature Science Foundation of China (Grants 41604011, 11673076, 11373001), the Geographic Information Engineering State Key Laboratory Open Fund of China (Grant SKLGIE2016-M-2-2), and the Outstanding Youth Fund of Information Engineering University (Grant 2016611705). The authors thank two reviewers for their comments on the manuscript and the editor for the handling of the review process.