Multi-parameter adjustment for high-precision azimuthal intersection positioning

The traditional azimuthal intersection method is viable for situations with only two control stations but a simple height averaging is not rigorous because the intersections vary in their distances from the two stations. In order to obtain the high-precision azimuthal intersections, this study presented a multi-parameter adjustment method, together with the Earth curvature correction and the atmospheric refraction correction models. This method is robust with varied distances between the control stations and the targeted intersections, without limitation of station quantity.• Based on the traditional space intersection, a multi-parameter adjustment model is added into the data processing for high-precision 3D positioning.• Both the Earth curvature error correction model and the atmospheric error correction model are included in the multi-parameter adjustment model, so the intersected points are more accurate than traditional intersections.

. a Sketch of the spatial azimuthal intersection method. A is one control Station, B is another control station and P is a targeted intersection such as a stake on the glacier. b Sketch of the earth curvature correction, the value of R is 6,378,137 m and the geodetic system is WGS84.

Azimuthal intersection method
In order to get the position of a target far away from the observatory, space intersection method is a good choice. All the azimuthal data can be collected by total stations or theodolites on the control stations. Using the classic space azimuthal intersection method, at least two control stations are needed to obtain the 3D coordinates of the target.
When two control stations, named A and B, are used to position one target, there are two pairs of horizontal and vertical azimuths from the stations at A and B. From these measurements, the 3D location of a stake can be calculated, and one redundant vertical azimuth exists, as shown in Fig. 1 a. Due to errors in the observations, the heights calculated from the vertical azimuths at both A and B differed, i.e., there was a vertical bias inside one intersection. Given that the locations of station A and station B were ( x A , y A , h A ) and( x B , y B , h B ), respectively, as acquired independently via GPS, and that the observed azimuthal records were ( α A , β A ) and ( α B , β B ) from stations A and B, respectively, the 3D location of a stake was traditionally calculated using Eq. (1) [5] .
In Eq. (1) , S A and S B are the horizontal distances from the intersection point P and the stations, The traditional space intersection method, as shown in Eq. (1) , is viable for situations with only two stations. Moreover, the S i varied from the stations. Therefore, a simple height averaging is not suitable for the error distributions at different distances, which means that the height averaging method is not rigorous.
In order to obtain high-precision location of the intersection point P( x P , y P , h P ), a multi-parameter adjustment model can be used to obtain the best solution. In this process, iterative calculations should be used, and the desired coordinates of the intersection point are set as the unknown parameters.
Based on the traditional space intersection, we used Eq. (1) to calculate the approximate intersected coordinates ( x, y, h ), and then to calculate the approximate horizontal and vertical azimuths via Eq. (2) .
Here, the horizontal azimuth has a scope of 0 ≤ α < 2 π , and the vertical azimuth has a scope of −π / 2 < β < π / 2 . Note that the value of α in Eq. (2) should be decided by coordinate y. When y − y i < 0 , α = 2 π − α i . Suppose the corrections of the coordinates of P are ( ˆ x , ˆ y , ˆ h ) , then the corrections of the azimuths, δ α i and δ β i , are as shown in Eq. (3) from a Taylor series expansion: Eq. (3) is the linear relationship between the coordinate corrections and the azimuthal corrections, which can be used for the programming, with consideration of the identical unit.

The parameter adjustment model
Since the errors in this study are mainly caused by the azimuthal observations, we can use the parameter adjustment model [2][3][4] to get the best solution, as shown in Eq. (4) : Here, the azimuthal correction vector In this study, the precisions of the horizontal and vertical observations were equal. Therefore, the weight matrix P was used as the unit matrix. The superscript T inside all the equations indicates the transposition of a vector or matrix. For example, the V T is the transposition of vector V , and the B T is the transposed matrix of B .

The residual vector of the azimuthal observations is
Here, the approximate values of the azimuths ( α A 0 , β A 0 , α B 0 , β B 0 ) can be calculated from Eq. (2) .
Note that the calculation in Eq. (4) requires iteration, which means the vectors B, L, and X should be revised every time. Consequently, the calculated values of the azimuths ( α A 0 , β A 0 , α B 0 , β B 0 ) should also be renewed with latest ( x, y, h ) after each iteration. Actually, when all the coordinate corrections are less than a predefined threshold value (e.g., 0.001 m), the iterations should be terminated, and the adjusted coordinates can be outputted.

Error corrections
The adjustment above reduces only the potential errors of the azimuthal observations. There were other errors that should be corrected, including the earth curvature correction, the atmospheric refraction correction and the zero-point reading correction. The earth curvature correction is especially remarkable in this research when the distances between the control stations and the intersections are over kilometres.
The earth curvature correction, δ R , used to amend the surface curvature (earth radius R ), causes errors in elevation, as shown in Fig. 1 b. The longer the plane distance S is, the greater the error is. The height error over a distance of 2.5 km reaches half a metre. The atmospheric refraction correction, δ air , is mainly used to reduce the height error caused by vertical atmospheric refraction, which may be several centimetres or even over ten centimetres in our research. Both δ R and δ air can be calculated according to Eq. (6) [6] .
After observing each target and collecting its azimuthal record, the station was supposed to be aimed back again to check the horizontal azimuth, creating a zero-point reading. Due to instrument errors and in situ environmental impacts, the zero-point readings are usually not zero. In this research, half the zero-point readings were used for the horizontal azimuthal observations while the other half was used for the intersection's 2D error evaluations. On the other hand, the zero-point readings can be entirely assigned to the observations or be evenly divided into the 360 °, but these two models are not recommended.
In this research, the coordinate reference systems is WGS84 and the value of R is 6,378,137 m. The vertical atmospheric refraction coefficient was fixed as k = 0 . 1 according to the results of previous works in the polar regions [6] . Since an error of 0.01 in k would lead to a bias of ~0.002 m in the elevation calculation, the potential error in k can be ignored.

Method validation
We have used this method to process the field azimuthal observations from Dalk Glacier, East Antarctica [1] . The two control stations were at two mountain peaks, and the targeted intersections were the stakes on the glacier. The stations' positions were acquired from GPS measurements, and the stakes' positions could be located with in situ azimuthal observations. After data processing using traditional space intersection method and our high-precision method, the comparison of position results showed that high-precision method corrected stake locations with a mean value of 0.156 m in height and 0.045 m in plane.