Application of Robust Estimation Method for Establishing 3D Combined Terrestrial and GNSS Network: A case of a Quarry in Lang Son, Vietnam

Recently, in Vietnam, the detection of geodetic measurements that contain rough errors as well as such data processing method has been considered as a key step in geodetic data processing, especially for large geodetic networks with many different types of measurements like 3D Global Navigation Satellite Systems (GNSS) network. On the other hand, mines in Vietnam often have complex terrains, so it is necessary to apply modern and flexible surveying methods in combination with ground and space measurements to build 3D coordinates control networks for management and exploitation to ensure sustainable development. Therefore, this research developed a Robust estimation method based on empirical weighting function for establishing 3D geodetic network combining terrestrial observation and GNSS vectors. The experiment on processing the combined network in Lang Son limestone quarry, Vietnam showed that the proposed method could be an effective solution for processing 3D terrestrial – GNSS geodetic network for mine surveying in Vietnam.


Introduction
Currently, as Vietnam is modernizing mining, the speed and volume of mineral exploitation are increasing in a wide range. Therefore, it is essential to evaluate the impacts of mineral exploitation affecting mine safety as well as to manage and forecast mineral deposits. In addition, the terrain of mine area is specific and complex, so it is necessary to study methods of building modern geodetic grids to meet multiple objectives such as monitoring of mine transfer, calculating exploitation volume, calculating mineral deposits, exploitation planning, etc.
The strong development of measurement technology with many modern high-precision measuring devices such as Global Navigation Satellite Systems (GNSS), total stations, etc. has contributed an effective solution in the construction of geodetic network, especially geodetic network for mining. However, each technique has its own drawbacks, in order to overcome the disadvantages of each measurement method, a common solution is to combine GNSS technology and traditional measurement technology to form a 3D terrestrial -GNSS geodetic network with high flexibility and accuracy (Krakiwsky and Thomson, 1974;Thomson 1976;Groten, 1977).
Geodetic data processing in general or adjustment of a geodetic network in particular is one of the important tasks in surveying. When constructing geodetic networks, measurement errors are inevitable. Therefore, network adjustment is required to find the most reliable values of unknown quantities. In this processing, rough error adversely affects the results of the adjustment problem. Modern adjustment theory has been studying the impact of rough errors on post-adjustment results and how to deal with them. Geodetic data obtained through statistics and analysis show that the prob-ability of appearing rough errors accounts for about 1% ÷ 10% (Tukey, 1962). One of the effective methods to handle rough errors is a robust estimation (Huber, 1964(Huber, , 1981Hampel, 1986;Koch, 2013Koch, , 2014Koch and Kargoll, 2013). Some scientists formulated weighted functions or various robust estimation methods such as the Danish method, the Tukey method, and the L1-norm method. However, the above methods consider the measurements to be independent and do not consider the correlation of the measurements in a geodetic network, especially the GNSS geodetic network.
Therefore, in this study, we use the Huber method and the weight function equivalent to two coefficients developed by Yang Y, et al. (2002) to develop a robust estimation with the application of empirical weight function to process 3D terrain -GNSS geodetic network data in order to meet the requirement of mine surveying in Vietnam.

Literature review of robust estimation 2.1 Huber method
The Huber method was proposed by (Hubber, 1981), and its function is expressed as follows: (1) where, v is the correction of the observation calculated from the previous iteration; k is constant, it is possible to take k=(2σ÷3σ), and the corresponding weighting factor is (2) Submission date: 06-03-2020 | Review date: 22-09-2020

Tukey method
According to Tukey (1962), Tukey function has the form as follows: Tukey's weighting function is as follows (5) where u=v/c.MAD, c is the regression coefficient.

Danish method
The principle of the Danish method is based on the indication of outliers by corresponding major corrections. After least squares adjustment, the priority weights of measurements are replaced by correction function. According to Kraup and Kubik (1983), the weighting function is as follows: (6) where c is constant and is often chosen as follows: c = 1.5 σ

L1 -norm method
The L1-norm is one of the most successful robust estimation (Harvey, 1993) and is expressed as follows (Wang et al., 2006): In L1-method, the corresponding weighting factor is (8) To solve the problem of weighting when v i =0, the weighting factor can be taken as (9) where k is a small number. The principle of correction is (10) On the other hand, with p i =p i w i , to be replaced in V = AX+ L to gain a standard equation system and its solution X is determined by the formula (11) as follows: (11) The above robust estimation methods consider the measurements to be independent and generally apply to small geodetic grids when the measurements have the same accuracy. However, for a large geodetic network with many types of measurement of different accuracy, especially 3D-GNSS geodetic networks, the correlation factor between measurements should be taken into account so that the above methods might be not effective

Robust estimation for 3D Terrestrial -GNSS network
To overcome the disadvantages of the Huber, Tukey, L1norm and Danish robust estimation method when applying for the 3D terrain -GNSS geodetic network, the robust estimation method for the measurements not of the same accuracy is used according to (Wang et al., 2006): Assuming that error equation is as following (12) where A nxn is a coefficient matrix, a 1 is the coefficient of matrix A, X tx1 is unknowns vector and needed to be determined, L nx1 is As measurements are not accurate, the matrix of weight P is as in the following form Replace V=AX+L in (17), then the standard set of equations of M estimation is (18) where: P is the equivalent weighting matrix, p i is the equivalent weighting element, w i is a weighting coefficient.
Parameter X of M robust estimation is determined as: According to Huber (1964), the general form of is ϕ where: c is the constant and normally takes the value of c=1.5÷2σ Independent measurements, the matrix of new weights is determined as: Dependent measurements, the matrix of new weights is determined as [14]: where γ ii , γ jj is the coefficient of weight loss . However, the method of using the two-factor equivalent weighting function [14] is only used in the GNSS grid. Therefore, to process the data of the 3D terrestrial -GNSS geodetic network including both ground and satellite observation, the proposed method is based on the Huber weighting function and the two-factor equivalent weighting function [14] to investigate and select the parameters of the function for independent measurements, dependent measurements, and c coefficient. Through experiments, the empirical weighting function is introduced for 3D terrestrial -GNSS geodetic network as follows: where c is constant and can be chosen as c = 1.5, σ vi is determined as the following formula: The Robust estimation includes four steps as follows: Step 1: n values of the geocentric coordinate system were calculated and transferred into n values of the geographic coordinate system ΔX, ΔY, ΔZ. The GNSS geodetic network is formed by edge vectors according to the principle of relative positioning. Measurements in the GNSS geodetic network are components of edge vectors ΔX, ΔY, ΔZ determined in the geocentric system (WGS-84) incorporated with their covariance matrix (WGS-84) symbolized as C: Simultaneously, the covariance matrix in the M system is determined based on the rotation matrix R and the covariance matrix in the geocentric system C XYZ : In which: Step 2: Adjustment of the network combining GNSS measurements (ΔX, ΔY, ΔZ) and terrestrial distance measurements S k =(k=1,2...n 2 ).
In this step, it is necessary to set up observation equations for GNSS edge vectors and terrestrial edges, including the coefficient matrix of each vector and free term L i : In which a T i is the coefficient of the correction equation of GNSS measurement value a S k ; is the coefficient system of correction equation of terrestrial edge measurements.
The general standard equation system is solved, an unknown vector X is gained as follows. (39) Step 3: Assessing the quality of the measurements after adjustment For accuracy assessment, root mean square error of weight unit σ is computed as follows: where n is the number of GNSS edges; n1 is the number of horizontal angles; n2 is the number of terrestrial edges; m is the number of unknown points.
Additionally, statistical standard χ2 is used to assess the reliability of the adjustment results according to the distribu-tion, if it meets the standard of χ2(χ2 r,1-α/2 ≤ VTPV≤ χ2 r,α/2 ) then adjustment results are accepted, if χ2 condition is not satisfied, the next step will be applied.
Step 4: If the observations contain coarse errors, the robust method with an empirical weight function computed by equation (26) is used to detect measurements that contain coarse errors.

Numerical example 4.1 Study area and materials
Lang Song province located at North East of Vietnam is a mountainous area, where there are many limestone quarries supplying raw materials to a large cement plant in Lang Son. The annual volume of limestone mining is very large; therefore, in order to serve the management of mining reserves, mapping the exploitation, setting up a precise geodetic network is essential. However, the topographical conditions of Lang Son quarry are complex (Figure 1), the placement of GNSS device at some control point is not satisfied due to unsatisfied threshold angle or near unsatisfactory power station. Therefore, to facilitate mine surveying, a 3D terrestrial -GNSS network was established. The combined network was illustrated in Figure 2. The network consists of 5 points, two original points are I, II, and three unknown points are III, IV, V. Due to the topographic conditions, GNSS device could not be placed at point I. Therefore, the edges I-III, I-IV, I-V were measured by an electronic total station (TS) with accuracy σ Si = (2mm+2ppm*Si) and edges II-III, II-IV, II-V, V-IV, II-IV were measured by GNSS technology with accuracy σ Si = (3mm+1ppm*Si). The data of TS edge measurement are presented in Table 1, the data of GNSS measurement in the geocentric coordinate system are shown in Table 2.

Experimental framework
The experiment includes five steps illustrated in below Figure 3. Where, Step #1. Reading GNSS and terrestrial data.
Step #2. Calculating and transforming measurement values of ΔX ΔY ΔZ and the covariance matrix from the geocentric coordinate system to the geographic coordinate system.
Step #3. Adjusting 3D terrestrial -GNSS geodetic network in case of the observations do not contain rough errors using the least-squares method. The adjustment results (i.e. correction vi) are considered as reference data for comparing to results of other Robust methods when the measurements contain rough errors.
Step #4. Assigning coarse errors: Firstly, coarse errors with the corresponding magnitude of 4σ, 5σ, 6σ,…, 9σ, 10σ (σ= 3mm) are randomly assigned for measurement of . Then the deviation of the corrections is calculated as Δv i =|v i -v i '| (39) where v i is the correction in case the measurements do not contain coarse errors; v i ' is the correction in case the measurements contain coarse errors.
Step #5. Coarse errors detection: The possibility of detecting coarse errors of experimental weighted function methods, L1-norm method, Tukey method, Danish method and Huber method is investigated. The conditions to detect observations containing coarse errors are presented in equation (26).

Results and discussion
In order to simplify the experiment, a computer module was programed using Microsoft Visual C language. This module can read observation data, both of GNSS and terrestrial data; estimate and transform variance and co-variance matrix, and adjust integrated network using LSM or Robust estimation method if the network considered as containing outliers. The main interface of this module is illustrated in Figure 4.
After adjusting the combined network without coarse errors using LSM, the coarse errors with different value were added into the measured value of ΔX v-iv . Then the robust estimation using the five methods was performed. The corresponding deviation of the obtained correction of ΔX v-iv comparing to the case without coarse errors was computed by equation (39) and shown in Tab. 3. Additionally, the deviations of Δv to the assigned coarse error were shown in Tab. 4.
It is clear that when applying robust estimation using experimental weighting function, the deviation of Δv to the assigned coarse error is smallest, and it gradually decreases when the assigned coarse error increases. Especially, it is closed to 0 when the assigned coarse error reaches 9σ (i.e. coarse error = 9x3mm = 27mm). It means that identifying coarse errors based on robust estimation using experimental weighting function is easier than the other methods.
In the second case, the assumption that 2 measurements are containing coarse errors, they are and , and the assigned coarse error value is 15mm. The corrections obtained by the five robust estimation methods were shown in Figure 5.
It can be seen in Figure 5 that experimental weighting method detected coarse error value approximates the assigned coarse error value (i.e. 15mm). This shows that, when using the experimental weighting function, the coarse error mainly affects the measurement value containing the coarse error and just slightly affects other measurements. Therefore, it is easy to identify the measurement containing coarse errors.
To assess the possibility of coarse errors detection by five methods, the percentage of observation containing coarse error detected was computed and Table 5. As can be seen in Table 5, the experimental weighting function method gives the best results compared to the L1-norm, Tukey, Danish and Huber methods, particularly when the measurement value contains errors less than 10% of the total number of measurements, the algorithm can detect 81-98%. The L1-norm, Tukey, Danish and Huber methods give lower results within a range from 50% to 86%. Therefore, when identifying rough errors for a 3D combined network, it is required to consider the correlation of GNSS measurements.

Conclusions
The experimental result showed that the Robust estimation method applying an empirical weighting function could quickly and accurately detect measurements containing coarse errors. Processing data of 3D terrestrial-GNSS network in mine area using experimental robust method can eliminate measurements that contain rough errors; as a result, it improves the accuracy of mine surveying such as determining mining reserves, monitoring mine slide. Therefore, it can be an effective solution for the processing of the integrated geodetic network data for mine surveying in Vietnam.
Establishing 3D terrestrial-GNSS network of Lang Son mine ease the network positioning in the national reference system; as a result, exploitation planning and management is more convenient. It is useful to establish 3D terrestrial-GNSS network.