Geodetic innovation in Chilean mining: The evolution from static to kinematic reference frame in seismic zones

2.1 History

In Chile, the georeferencing of mining concessions in the cadaster has traditionally relied on two classical geodetic reference frames:

• South American Provisional Datum 1956 (PSAD56) from the north of the country to the parallel –43°30′.

• South American Datum 1969 (SAD69) from parallel –43°30′ to the south of the country.

Before the adoption of these cRFs, the georeferencing of mining concessions was carried out using arbitrary local coordinates of astronomical origin, which ceased to be utilised due to their local and limited nature. The geodetic infrastructure for PSAD56 and SAD69 was established between the 1950s and 1970s by the IGM, the agency responsible for generating the national cartography (Ministerio de Defensa Nacional, 1963). The IGM and Inter-American Geodetic Survey (IAGS) conducted observation campaigns to materialise PSAD56 and SAD69 by establishing national geodetic networks (Instituto Geográfico Militar, 2020). The campaign finished in mid-1970, but due to a series of large-magnitude earthquakes since 1960, including the Valdivia mega-earthquake in 1960 (9.5 Mw; Astroza and Lazo, 2010), it was necessary to make re-observations in the areas affected by the earthquake (Instituto Geográfico Militar, 1970). The static definition of PSAD56/SAD69 and partial updates of the network (Instituto Geográfico Militar de Chile, 2005) have generated two deteriorated cRFs that were quickly outdated and had low applicability due to the geodynamic conditions of Chile. Therefore, in 2002, the IGM started transitioning from these cRFs to mRFs such as SIRGAS. The coordinate couples between the cRF and the mRF enabled a connection to be generated between them from transformation parameters. As preliminary background, in Chile, there is only the IGM project, which used the transformation parameters calculated by the National Imagery and Mapping Agency (NIMA) for the transition from PSAD56/SAD69 to WGS 84 (World Geodetic System 1984) with European Petroleum Survey Group (EPSG) codes 6971, 6972, 6973, and 6977. WGS 84 is equivalent to SIRGAS at a centimetric level for the observation epoch. With this, the IGM suggested that SIRGAS could serve as Chile’s official reference system in 2003. IGM published the transformation parameters between PSAD56/SAD69 and SIRGAS to adopt an mRF (static) in Chile; these parameters are only for cartography, not mining (Instituto Geográfico Militar de Chile, 2008). Due to the IGM’s proposal, the SERNAGEOMIN identified and analysed different methodologies to implement the new reference frame and to transform all existing information from PSAD56/SAD69 into SIRGAS; however, it did not carry out any of them due to a lack of tools or viable solutions (Banderas Briceño, 2010). Our scientific proposal aims to fix the previous problem.

2.2 Data

Table 2 shows the data used for this research:

The cRF coordinates are obtained directly from the historical database. In contrast, the REDGEOMIN coordinates are obtained from scientific processing from the RINEX files of stations with open data available in Chile. The processing will be seen in the following section.

3.1 REDGEOMIN processing strategy

REDGEOMIN is the first Global Navigation Satellite System (GNSS) geodetic network designed for mining purposes that integrates all agencies with active and open access geodetic infrastructure in Chile, as suggested by the United Nations (UN-GGIM, 2015). The network (Figure 4) consists of approximately 390 Continuously Operating Reference Stations (CORS) distributed throughout the country and 30 IGS stations that aligned the mRF to IGb14 (Figure 5) restricting to the corresponding IGS weekly positions (Brunini et al., 2012).

Figure 4

REDGEOMIN Stations. Figure adapted from IAG Beijing (Tarrío Mosquera, Inzunza, et al., 2021). The figure was created using GMT software (Wessel et al., 2019).

Figure 5

IGS stations used to fix REDGEOMIN. Adapted from IAG Beijing (Tarrío Mosquera et al., 2021). The figure was created using GMT software (Wessel et al., 2019).

The network calculation was performed using two scientific software, BSW (Bernese Software) (Dach et al., 2015) and GG (Gamit GlobK) (Herring et al., 2015), and processed according to IGS (Johnston et al., 2017) and SIRGAS (Tarrío Mosquera et al., 2021) guidelines. BSW’s processing is the same as the SIRGAS-CON analysis centres (Sánchez et al., 2022). It makes use of 7-day observations (Sunday to Saturday) and metadata such as precise ephemerides, antenna files (ANTEX), rotation parameters, tropospheric delay and ionospheric models, and tide files. During processing, the phase ambiguities solution is performed according to the baseline length: (1) direct solution using L1 and L2 for baselines from 0 to 20 km, (2) using L3 and L5 for baselines from 20 to 200 km, (3) wide lane for baselines from 200 to 9,000 km, and (4) near ionosphere-free for baselines from 18 to 5,600 km normal standard (Dach et al., 2015). Daily normal equations are calculated using the double difference method. Baselines are constructed based on the maximum common observations of associated stations. The normal equations are obtained day by day once the ambiguities are resolved. The normal equations are combined with daily loosely constrained solutions on the central day (Wednesday), resulting in a weekly solution (Brockmann and Gurtner, 1996). Finally, the weekly solutions are aligned to IGb14 with IGS weekly positions, considering 30 IGS stations as fiducial points. The selection of these stations was made considering (1) their belonging to IGS/SIRGAS, (2) the time series of the station with high stability and low variability, and (3) their proximity to Chile. IGS/SIRGAS stations outside Chile are used to align the solutions to the IGS reference frame, while those within Chile are used to monitor existing deformations. This selection method is applied in processing regional networks, such as EUREF (Altamimi 2003; Altamimi et al., 2004; Boucher and Altamimi, 1992). The processing with GG was similar.

Four annual solutions were obtained from 2019 to 2022 to analyse the RF annual kinematic. The objective is to have transformation parameters that relate PSAD56/SAD69 to REDGEOMIN in the project’s most immediate period to the current one, 2022, in this case.

There are CORS before 2019 in Chile with intermittent data distribution since 2009 and heterogeneous distribution. These data were used to calculate the deformation model and the coordinates of the CORS that materialise REDGEOMIN. ADELA comprises a time series generated from data from 2009 to 2022, and it includes a thin plate spline (TPS) interpolation algorithm for areas with no CORS. The coordinates of the time series are obtained from REDGEOMIN weekly processing and supplemented in regions without data with precise point positioning (PPP) coordinates generated by the Nevada Geodetic Laboratory (NGL) (Blewitt et al., 2013). The time series allows the generation of a composite function of the trajectory of each station. The function process is modelled with a joint model considering interseismic, coseismic, and postseismic periods (Bevis and Brown, 2014). The postseismic deformation (PSD) is modelled with logarithmic and/or exponential functions (Sobrero et al., 2020). The whole process is performed according to equation (1):

where

(1) ∑ i = 1 n p + 1 pi ( t − tr ) i − 1 is the linear velocity,

n p is the order of maximum power of the polynomial, t is the time of each epoch, and tr is the reference epoch.

(2) ∑ j = 1 n j b j H ( t − t j ) , earthquakes, and equipment changes,

n j is the number of jumps (Heaviside), b j is the j th jumps, and t j is the epoch at the jump.

(3) ∑ k = 1 n F [ p a , k cos ( ω k ( t − t 0 ) ) + p b , k sin ( ω k ( t − t 0 ) ) ] , periodic signal,

n F is the number of frequencies used to model the annual displacement cycle, and ω k , 2 π / τ k , terminology associated with cycle periods:

(4) ∑ i = 1 n log a i log 1 + ∆ t T i + ∑ i = 1 n exp b i [ 1 − e − ∆ t T i ] , PSD

a i is the amplitude of the transient, T i is the time scale, and ∆ t is the time since the earthquake occurred.

Incorporating PSD values is necessary since the linear and/or periodic component fails to model the geodynamics of Chile (Sánchez and Drewes, 2020; Gomez et al., 2015). The complexity of crustal movements in Chile is shown by the CORS displacements (red arrows). Figure 6 (left) shows the crustal behaviour in the interseismic and coseismic periods, and Figure 6 (right) represents the postseismic displacement. Both figures display the deformation heterogeneity throughout the country, which complicates using a modern static RF (Tarrío Mosquera et al., 2021; Blick et al., 2014).

Figure 6

Displacement in interseismic and coseismic for Iquique 2014 (8.2 Mw) and Coquimbo 2015 (8.3 Mw) earthquakes. Postseismic displacements (2019–2021). Figure adapted from IAG Beijing (Tarrío Mosquera et al., 2021). The figure was created using GMT software (Wessel et al., 2019).

In areas without active stations, the benchmark displacement needed to be interpolated from the information obtained from the time series. For this purpose, the model uses the TPS interpolation algorithm to interpolate the coordinates (Keller and Borkowski, 2019). This interpolation allows us to obtain the geodetic displacement of the station directly without modelling the movement as it would be, for example, with a least-squares collocation. Figure 7 exemplifies the latitude, longitude, and height motion due to the devastating earthquakes on Illapel 2015 (8.3 Mw) and Iquique 2014 (8.2 Mw). Figure 8 shows the coseismic displacement for the above earthquakes.

Figure 7

TPS interpolation for the earthquakes. IAG Beijing (Tarrío Mosquera et al., 2021).

Figure 8

Mapping coseismic deformation. The figure was created using GMT software (Wessel et al., 2019).

Figures 7 and 8 show metric deformations with divergent directions for nearby CORS (approximately 50 km). The exposed methodology allows obtaining an mRF, including the seismic variable, and modelling the divergence in areas without stations. The results are an mRF with (1) REDGEOMIN coordinates from epoch 2019.00 to 2022.00, (2) time series 2009–2022, and (3) the ADELA deformation model. An online calculator (in version beta) has been developed using this information, available at https://adela.usach.cl/.

3.2 Passive network adjustment

The PSAD56/SAD69 coordinates for the passive geodetic markers forming the secondary network were obtained from different sources, mainly from 1965 to 1970 (Instituto Geográfico Militar, 1965, 1970). GNSS measurements permit to compute the coordinates concerning REDGEOMIN of 158 benchmarks of the native PSAD56/SAD69 network of the IGM (Figure 9). In addition, observations for 31 benchmarks were also provided by SERNAGEOMIN. The epoch of these coordinates, before 2019, was updated to 2022 with the ADELA deformation model.

Figure 9

PSAD56/SAD69 measured points. Figure adapted from REFAG 2022 (Tarrío et al., 2022). The figure was created using GMT software (Wessel et al., 2019).

This passive network was measured in three phases corresponding to the years 2019, 2020, and 2021–2022. The observations were processed at epoch 2022.00 to obtain updated coordinates at the epoch closest to the project implementation. The adjustment was made to tie the secondary network observations to CORS from the primary network. The classic network’s points have a decimetric precision. Therefore, in order not to worsen the RMS transformation, the passive points’ precision was established in centimetres. The secondary network was processed using commercial and non-scientific software due to the type of observation and the reasons mentioned above (Hamidi et al., 2017; Schütz, 2017). Table 3 shows the processing parameters used for the secondary network.

3.3 Datum transformations

The PSAD56/SAD69 coordinates were revised before calculating the transformation parameters. The outliers were evidenced using a simple difference of geodetic coordinates in both systems (PSAD56/SAD69-REDGEOMIN). Figure 10 shows the 11 points with PSAD56/SAD69 differences greater than 3σ (red arrows), whose coordinates were eliminated for the transformation.

Figure 10

PSAD56/SAD69 outliers. The figure was created using GMT software (Wessel et al., 2019).

In this section, there is a limitation by law, which implies that the concessions must not be deformed in the new system when transformed from the old one, which means that, technically, a 2D transformation (T2D) only with translations must be used. The T2D will be performed per the UTM zone (Table 1). However, such transformation has significant limitations in a region with no linear deformation as it is the entire territory of Chile. Consequently, it was decided to evaluate the use of another transformation method that could permit a better adjustment between both systems, namely a three-dimensional Helmert transformation (seven parameters) (H7P) and two NTv2 grids. One of the advantages of the NTv2 approach over classical conformal transformations (T2D and H3D) is that it makes it possible to employ subgrids in more bounded areas to add to the primary grid (Garnero, 2014) without affecting the accuracy of the surrounding concessions. The number of transformations evaluated was determined based on the residue values to analyse the deformation of the cRF and the number of points exceeding the residual threshold.

4.1 REDGEOMIN assessment

The evaluation of REDGEOMIN results is performed by analysing the internal and external consistency as follows (Brunini et al., 2012; Costa et al., 2012):

1. The mean standard deviation for station positions after aligning the network to the IGS RF indicates the formal error of the final combination. The weekly coordinate repeatability after combining the individual daily solutions provides information about the internal consistency of the combined network.

2. Comparison with the IGS/SIRGAS weekly coordinates indicates external consistency with the IGS global network and SIRGAS-CON.

Table 4 presents the weekly repeatability of the four computed REDGEOMIN annual solutions. The values for internal consistency are less than 2.3, 2.3, and 5.3 mm in E, N, and U, respectively (in topocentric systems or local systems), with a confidence level of 1σ (Dach et al., 2015). The E, N accuracy remains at about 2 mm for all epochs, while U decreases from 2019.00 to 2022.00 from 3.6 to 5.1 mm. This increase may be because from 2019 to 2022, the REDGEOMIN stations increased by approximately 30 CORS.

Figure 11 shows that the solutions are aligned to IGS with an accuracy below 1 mm in the three components for all campaigns. In contrast, Figure 12 indicates that the solutions are aligned within SIRGAS with 1, 1, and 2.5 mm accuracy in E, N, and U, respectively. The accuracy improvement was because, in epoch 2022, there were more IGS/SIRGAS stations with data to process.

Figure 11

IGS-REDGEOMIN alignment.

Figure 12

SIRGAS-REDGEOMIN alignment.

The current cartographic RF for Chile (SIRGASChile) is available in epoch 2021.00, and it has accurate transformation parameters for mapping, specifically at the 1:25,000 scale (Rozas et al., 2021). The alignment results are shown below to evaluate the interoperability of REDGEOMIN and SIRGASChile. Figure 13 shows the IGS/SIRGAS/SIRGASChile alignment versus REDGEOMIN 2021.00, where it is observed that the largest deviation is presented when calculating the alignment between REDGEOMIN and SIRGASChile.

Figure 13

IGS-SIRGAS-SIRGASChile v/s REDGEOMIN.

There is not much information about SIRGASChile’s processing strategy; according to Rozas et al. (2021), we believe that this worsening of accuracy could be due to fewer IGS and SIRGAS fiducial points and to different scientific processing models employed. In verbal conversations with IGM, colleagues indicate that they use a Bernese reduced script for processing the Chilean network (SIRGASChile) (Parra et al., 2014). The differences may, therefore, be due to including fewer fiducial points and a different processing strategy.

ADELA was calculated from REDGEOMIN weekly solutions; for this reason, and to evaluate the accuracy of these, we performed a cross-validation of the modelled positions. Figure 14 represents the 22 stations used for validation; the accuracy is σ = ±5 mm, with a confidence level of 1σ (dashed red line), and Figure 15 shows the location of these stations along Chile.

Figure 14

CORS statistics used for ADELA cross-validation.

Figure 15

CORS distribution used for ADELA cross-validation. The figure was created using GMT software (Wessel et al., 2019).

Therefore, the model’s accuracy is according to mRF.

4.2 Passive network accuracy

The passive network consists of 189 points, of which 11 were eliminated, leaving 178 points, as shown in Figure 16. Table 5 shows the RMS values resulting from the processing and adjustment of the secondary network. The horizontal and vertical components’ accuracy is less than 1 and 2.5 cm, respectively. The distribution of these errors (frequency vs residuals) is also shown in Figures 17 and 18.

Figure 16

Passive network. The figure was created using GMT software (Wessel et al., 2019).

Figure 17

Precision EN vector processing.

Figure 18

Precision Up vector processing.

Figures 17 and 18 illustrate the precisions of the points are distributed within a range. The precision of the processed points (cm) is below that of the classical system points (dm). Therefore, the precision of the transformation will also be in the decimetre order in the best case.

4.3 Transformation parameters

The transformations used were various: conformal and transformation grids. In the first case, a T2D that complied with the mining law was also used but had an RMS almost triple that of the others because it does not absorb the deformations of a seismic country like Chile.

4.3.1 T2D transformation

About 158 of the 178 secondary network points have PSAD56 coordinates, while the remaining 20 have SAD69 coordinates. The transformations were calculated by zone, aligning with the law requirements (SERNAGEOMIN, 2019). This decision stems from the fact that the coordinates of mining concessions are recorded in the UTM system, as mandated by the mining code (Ministerio de Minería, 2014). Table 6 shows the RMS values for the T2D transformation using different thresholds for the maximum value of the accepted residuals. For PSAD56 UTM zones 18 and 19, four T2D transformations were performed. The maximum value of the accepted residuals varies from 1.5 to 3.6 m in zone 18 and from 1.0 to 3.0 m in zone 19. For SAD69 UTM zones 18 and 19, two T2D transformations were performed. In this region, the maximum value of the accepted residuals varied from 2.5 to 5.0 m in zone 18 and from 1.0 to 2.0 m in zone 19. These different values were selected considering the number of existing control points and the values of the estimated residuals.

Table 6

Translation 2D transformation parameters

PSAD56-REDGEOMIN T2D									SAD69-REDGEOMIN T2D
Zone 19 S									Zone 19 S
Maximum residuals	3.00 m		2.00 m		1.50 m		1.00 m		Maximum residuals	2.00 m		1.00 m
Total points	137								Total points	15
Used points	116		97		72		46		Used points	13		11
∆E (m)	−184.43	±1.04	−184.41	±0.94	−184.15	±0.66	−184.04	±0.47	∆E (m)	−68.69	±0.79	−68.41	±0.43
∆N (m)	−374.47	±1.25	−374.63	±0.95	−374.65	±0.79	−374.56	±0.59	∆N (m)	−19.26	±0.29	−19.23	±0.31
RMS (m)	1.15		0.94		0.72		0.52		RMS (m)	0.57		0.35
Zone 18 S									Zone 18 S
Maximum residuals	3.60 M		3.00 M		2.00 M		1.50 M		Maximum residuals	5.00 M		2.50 M
Total points	23								Total points	5
Used points	16		11		10		7		Used points	4		3
∆E (m)	−237.90	±1.90	−237.27	±1.21	−237.08	±1.08	−237.52	±0.96	∆E (m)	−86.67	±3.74	−88.32	±2.15
∆N (m)	−362.00	±2.21	−360.78	±1.39	−360.51	±1.11	−359.98	±0.76	∆N (m)	−26.79	±2.14	−26.15	±2.09
RMS (m)	1.99		1.24		1.04		0.80		RMS (m)	2.64		1.50

Table 6 and Figure 19 show the point distribution in the study area for PSAD56-REDGEOMIN zone 19. The uncertainty improves with a lower threshold value for the maximum accepted residuals, as would be expected. However, this means that larger groups of points are excluded. In fact, there are only a few points in the southern part, implying that the estimated transformation parameters have greater uncertainty in that area. According to Table 6, the optimal parameters are those with the most points because the value of the parameters is practically the same, but more representative points are used.

Figure 19

T2D PSAD56-REDGEOMIN-zone 19 distribution. The figure was created using GMT software (Wessel et al., 2019).

We also compared our computed parameters with the estimated by NIMA (IGM adopted for the transition from PSAD56/SAD69 to SIRGAS). NIMA transformations are divided by zones (from one latitude to another); therefore, to transform PSAD56 coordinates, the parameters that go from 17°30′S to 26°S, from 26°S to 36°S, and from 36°S to 43°30′S were used. In the case of SAD69, NIMA only has transformation parameters that cover the area of Tierra del Fuego (south of Chile). Table 7 shows the statistics obtained by applying the NIMA and T2D transformations. Table 6 shows that most transformations’ average and standard deviations are lower when the transformation parameters calculated in these solutions are applied than when they are obtained using the NIMA parameters.

Table 7

NIMA vs REDGEOMIN parameter comparison

Residuals	E (m)	N (m)	E (m)	N (m)
	NIMA/IGM EPSG:6971-6972-6973/10135-10136-10137		T2D<3.0 m PSAD56-REDGEOMIN zone 19 S
Max	1.81	5.30	2.11	3.22
Min	−4.19	−5.68	−2.20	−5.93
RMS	1.18	2.21	1.02	1.97
Average	−0.55	−0.49	−0.02	−0.56
	NIMA/IGM EPSG:6972-6973/10136-10137		T2D<3.6 m PSAD56-REDGEOMIN zone 18 S
Max	6.57	6.42	2.58	3.03
Min	−1.36	−3.96	−4.99	−6.58
RMS	2.36	3.23	2.06	3.23
Average	2.15	0.95	−0.04	−1.70
	NIMA/IGM EPSG:6977/10138		T2D<2.0 m SAD69-REDGEOMIN zone 19 S
Max	1.43	0.76	1.27	0.64
Min	−1.35	−0.22	−1.68	−0.58
RMS	0.81	0.27	0.79	0.29
Average	0.28	0.30	0.00	0.00
	NIMA/IGM EPSG:6977/10138		T2D < 5.0 m SAD69-REDGEOMIN zone 18 S
Max	−8.53	−9.22	4.95	2.79
Min	−17.56	−13.52	−4.05	−1.92
RMS	3.74	1.98	3.74	2.14
Average	−13.47	−11.81	0.00	0.00

4.4 Additional solutions

As previously mentioned, other transformations were also computed to analyse the cRF deformation. Tables 8 and 9 statistics result from the H7P PSAD56/SAD69-REDGEOMIN tested transformations; these are calculated indistinctly of the zone. The threshold for the maximum accepted residuals ranged from 1.5 to 4.0 m.

Figures 20 and 21 show the distribution of accepted points in the study area for PSAD56-REDGEOMIN and SAD69-REDGEOMIN, respectively. The standard deviations, especially in SAD69, are high, possibly due to the distance to the datum point (Chuá, Brazil). Another aspect may be that the H7P uses ellipsoidal height in cRF, but we only had an orthometric height. It is important to mention that there are no accurate geoid undulations available in cRF, so we had to iterate the ellipsoidal height to obtain the orthometric one (González Matesanz, 2012).

Figure 20

H7P PSAD56-SIRGAS/REDGEOMIN distribution. The figure was created using GMT software (Wessel et al., 2019).

Figure 21

H7P SAD69-SIRGAS/REDGEOMIN distribution. The figure was created using GMT software (Wessel et al., 2019).

As with the T2D transformations, the H7P transformations also do not properly fit much of the study area when the threshold value for accepted residuals decreases.

Chile is a long (4,000 km) and narrow (200 km average) country, with an enormous deformation that the cRF does not show. The T2D and H7P transformations cannot absorb that singularity, which is unique worldwide. Therefore, NTv2 grids were evaluated. Table 10 shows the statistics for calculating the PSAD56 and SAD69 NTv2 grids. In both cases, the threshold value for the maximum value of the residuals was 1.0 m.

A scientific proposal must have applicability, so we tested the NTv2 in the mining cadaster. Therefore, we assess the impact of the change in the surface of the concessions using the NTv2 instead of T2D. Table 11 and Figure 22 show that the areas of Atacama, Libertador Bernardo O’Higgins, Maule, Ñuble, Bio Bio, Araucanía, and Los Lagos present more significant variations in the size of the concessions.

Figure 22

Areas with the highest percentage of variation in concessions when applying NTv2 grids. The figure was created using GMT software (Wessel et al., 2019).

The area throughout Chile shows no significant variation, allowing for transformation to REDGEOMIN@2022.00 with more accuracy than T2D.