Optimization of deformation monitoring networks using finite element strain analysis

1 Introduction

The determination of changes in the shape of man-made constructions or the surface of the earth is of importance in geodesy, where we need to carefully monitor deformable objects. The purpose of deformation monitoring is to study the behavior of a deformable object in short or long time intervals and estimate the possible displacement or deformation within the object. Basically, these studies are performed to prevent disasters due to earthquakes and landslides as well as the progressive destruction of large-scale structures such as dams and bridges. Surveying engineering provides discrete information such as displacements of the network points, while deformation is a continuous phenomenon. Therefore, an approach should be developed to estimate the displacement field all over the object.

It is quite common in some of the previous studies about deformation analysis to consider the whole deformable body as one object and compute the deformation parameters, i. e. strain and differential rotations, of that (see for instance: Doma [9], and Kuang [13], pp. 292–300). Thus, the obtained parameters describe the deforming condition of the whole body. However, a simple solution to find a more precise deformation model through a deformable object is to split the body into small pieces, i. e. finite elements [7]. The procedure can be performed with the help of the Delaunay triangulation [14], where the body is split into small, non-overlapping and possibly equilateral triangles with the vertices of known points. The advantage of this approach is to check the quality of the least squares fitting of the displacement model to deformation parameters of each element, and consequently perform a meticulous investigation of the displacement field in each element. For instance, Welsch [21] studied this method to determine the horizontal strain patterns from geodetic observations, and Kiamehr & Sjöberg [12] investigated the three-dimensional finite element method to analyze the surface deformation patterns of the Skåne area in Sweden. Finite Element Strain Analysis (FESA) is an approach to analyze the deformation of a deformable body and could be either performed on triangles and quadrangles in the plane or in the three-dimension (3D) by considering tetrahedral and hexahedral elements.

The precision of deformation parameters that are obtained by the FESA approach is used here to define an appropriate precision criterion matrix. This matrix acts as the core of the optimization procedure, where the goal is to fit the current Variance-Covariance (VC) matrix of the network to the criterion. Optimization of a deformation monitoring network helps us to design an observation plan, which fulfills all the pre-defined and required criteria for the network such as precision.

In this paper, we studied a Global Positioning System (GPS) monitoring network in Skåne in southern Sweden and developed a methodology to optimize that network. This network has already been used in another study to investigate the effect of GPS baseline correlations in designing an optimal observation plan for the network [2]. However, in this study, the main problem is how to involve the precision of deformation parameters of each finite element in design and optimization of the Skåne GPS network. Although the activity of this network was only limited to the years between 1989 and 1998 and it is not active anymore, this study is conducted to motivate future measurements with an optimal observation plan, which would enable the use of the existing network to detect possible deformations with a fewer number of static GPS observations.

2 Deformation parameters

Deformation can occur in an object due to different external forces and cause changes in shape, dimension, and position of the object. The role of time is imperative to be taken into account in deformation subjects. Most of the changes in the appearance of an object are found as time elapsed. Therefore, in surveying a deformable body for monitoring purposes the measurements are repeated epoch-wise or continuously in time.

If an object undergoes some homogeneous deformations, the parallel and straight lines of the object will remain unchanged, while in non-homogeneous type, those lines will become non-parallel and curved. To simplify the deformation model, it is very common to consider a homogeneous deformation and strain all over the region. It should be noted that in case of small-scale objects, it is reasonable to deal with such assumption, where almost the whole body has a similar behavior against the exerted force. On the other hand, in a large-scale deformable object, (e. g. continents) it is required to consider non-homogeneity of the ground in establishing a geodetic network. It means that a single deformation model cannot fit the displacement field of the whole area and explicate the occurred deformation.

The concept of strain in a geodetic network can be mathematically defined as the rate of change in the displacement field (ξx,y,z=α,β,γT) of an object with respect to position (x=x,y,zT), where α, β and γ represent displacement components in the directions of x, y and z, respectively. The strain matrix S, which consists of linear displacement gradients, is written as [20]:

Decomposition of S into symmetric and skew-symmetric matrices yields the translation T and rotation Ω matrices, respectively. We can write [3]:

Applying the derivatives in Eq. (1) in (2) yields:

where the symmetric strain tensor T consists of the normal (εx, εy, εz) and shear (εxy, εyz, εxz) strains at a point, and matrix Ω describes the deformation rotations (ωx, ωy, ωz) around three axes at a point. All these parameters, i. e. normal and shear strains as well as deformation rotations are known as deformation parameters. The normal strain scalars express the expansion and contraction. It should be noted that the strains describe only the relative displacement of points, and have not been affected by translation of the deformable object. The principle strains are the eigenvalues of the strain tensor and the directions of principle axes are the eigenvectors. Negative principles depict contraction versus expansion by positive principles [20].

Briefly, any homogeneous deformation has three components: strain, rotation, and translation. Generally, these components have non-zero values, but not necessarily. An object is said to have pure deformation if only one of these components becomes non-zero. For instance, an object has a pure rotation, if there are no strain and translation in that block (Means [15], p. 145).

3 Deformation models

Deformation surveys provide a priori information from different measuring techniques to build a proper deformation model for an area and analyze the corresponding deformation components. A number of typical deformation models in two-dimensional space are introduced and discussed in Chrzanowski & Secord [5] and Setan & Singh [19]. Chrzanowski, et al. [4] developed a generalized approach based on the least squares method to fit a deformation model to the displacement field obtained from repeated measurements of deformations. To estimate the deformation parameters, the whole area is generally considered as a non-continuous deformable body consisting of separate continuous deformable blocks. Due to the vastness of the study area and its geographical location, the forming blocks may undergo relative rigid body displacements and rotations, so each block has its own deformation model (Kuang [13], pp. 175–191).

For monitoring deformations by geodetic methods, usually, a network is established to cover the deformable body. The networks can be built either as a relative or a reference network according to the purpose of monitoring. In the former type, all the network points are located on the deformable body, while in the latter, it is assumed to set up some points outside of the deformable body as a reference point to determine the absolute displacement of the object points [5].

Deformation equations in a 3D space can be written for any point i located within the deformable body as:

and

In order to separate the nine unknown deformation parameters (strain and differential rotation components) from known positions (x, y and z) in a model, Eqs (4) and (5) can be reformulated in a matrix form as:

The shortened form of Eq. (6) is written as:

where d is a displacement vector, B is a coefficient matrix of deformation parameters, and e represents a vector which contains the deformation parameters. This equation shows that the displacement field can be approximated by fitting a selected deformation model to displacements determined at discrete points. Actually, fitting a plane surface to each displacement field, results in a very simple determination of the strain components.

4 Finite element strain analysis (FESA)

Once we want to estimate and analyze the deformation parameters in an area, we need to consider the nature of geodetic networks, which is discrete. The only available information from these networks for deformation estimation is the displacement vector of the network points. The piecewise polynomials presented in Eq. (5) can be used to interpolate the discrete point displacements and provide continuous information for a deformable body [8].

Additionally, assuming a homogeneous strain field, where the strain is evenly distributed throughout a deformable body is often not desirable in large-scale objects and regions [21]. Therefore, to get an insight into the local variety of strain patterns, the area of investigation can be fragmented into smaller elements where deformation components are constant. The simplest choice for the subareas is triangles or tetrahedrons of geodetic benchmarks with no overlay and gap between adjacent elements. The best triangulation is the one with optimal angles and as near equilateral as possible, i. e. Delaunay triangulation, which provides the largest angle vector and optimizes the connectivity for a given set of points [18]. The FESA enables us to firstly, obtain deformation information for each individual element and secondly, use the attained a priori information as a basis for defining a criterion matrix for the optimization procedure.

Grafarend [11] showed that the tetrahedron is a fundamental geodetic figure in order to determine the characteristic coefficients of the deformation equations, i. e. strain and differential rotation components in a 3D space.

After splitting the deformable body into finite elements by the 3D Delaunay triangulation technique, one may associate the illustrated tetrahedron in Fig. 1 with any of those elements. The vertices of the tetrahedron (a, b, c and d) represent any four irregularly located points in the network, which form an element with respect to the Delaunay optimal conditions.

Figure 1

A tetrahedron as a 3D finite element.

In order to be able to determine nine unknown deformation components in each element by the least squares, we need to write three deformation equations for all four points of the element. Considering Eq. (7), estimation of the deformation parameters is based on the following adjustment model:

where Δl is an observation vector containing the differences of coordinate measurements in two successive epochs with the residual vector of v. The element design matrix D describes how the four points in each 3D element are connected with six lines. Hence, the dimension of this matrix for each tetrahedral element j will be 6×3 by 4×3. B is the coefficient matrix of the deformation model, and e is the vector of unknown deformation parameters.

In case of dealing with more than one element, the matrices D and B may contain the individual coefficient matrices from all elements. Therefore, the structure of these matrices can be shown as:

where

and

where

The blkdiag ∙ operator makes a block diagonal matrix. n is the total number of fragmented elements, and j represents the number of each element.

Finally, all unknown deformation parameters of finite elements (ej) are stacked below each other in the unknown vector e. There are totally nine unknown parameters in each ej: six strains and three differential rotations

where

The estimated unknown vector of deformation parameters eˆ from Eq. (8) and its corresponding VC matrix Ce assuming two epochs of observations are as follows (Kuang [13], p. 262)

and

5 Optimization procedure

The main purpose of developing the optimization procedure for a geodetic network is to design a precise and reliable network that excludes all unnecessary observations from the observation plan. The principles of the network design was introduced by Grafarend [10] as the zero-order design to seek the best network datum, the first-order design to choose the best location for the network points, the second-order design to find the optimal weights for the observations, and the third-order design to add extra points or observations to the network. In this paper, we deal with the second-order design problem to remove unnecessary baselines from the observation plan, so the network can fulfill the demanded quality criteria (precision and reliability) even without these baselines. Basically, the deformation parameters in relative monitoring networks are datum independent as the point displacements are obtained relative to each other at different observation epochs [6]. Therefore, the concept of the zero-order design is of no practical importance in such networks. The first-order step is also not performed in this study due to minor effects on network configuration in providing the quality criteria of a Global Navigation Satellite Systems (GNSS) network.

In order to acquire the requested precision of the network in the second-order design stage, we need to properly define a criterion matrix. In this study, we define a criterion matrix for displacement of network points (Cdc) based on the precision of deformation parameters. Recalling Eq. (7) and applying the error propagation law to this equation, the covariance matrix of the displacements can be obtained as:

The Bc matrix in Eq. (17) has the same structure as the B matrix in Eq. (11), if we write the deformation equations separately for each element. However, in the optimization procedure, we need to tie up the elements with the geodetic observations. Therefore, the Bc matrix will have a different structure for the network points that are involved in more than one tetrahedron element. In other words, to hold the consistency in dimensions of the VC matrix of the GNSS observations and the criterion matrix, the average of the deformation equations are used for points with more than one element. The configuration matrix of deformation parameters Bc consists of a number of block matrices Bibc as

where xi, yi and zi are the coordinates of network point i, and n′ is the number of elements related to point i. For instance, if we assume that point number 1 is used to generate 3 elements out of totally 5 existing elements, i. e. element number 2, 3 and 5, then the structure of B1c will be as follows

To form the matrix B1c in Eq. (19), the three deformation equations (α1, β1 and γ1) are written for the point number 1 in each of these five elements. It is obvious that there are nine deformation parameters in each element. Thus, the dimension of this matrix would be 3 by 9×5. Eventually, the matrix Bc can be structured as

where m is the number of network points, and n is the total number of elements.

Moreover, the existing VC matrix of the displacements in a geodetic network can be expressed as (Kuang [13], p. 262):

where σ02 is the a priori variance factor that is predominantly assumed to be 1 in the design stage, A and P are respectively the configuration and weight matrices of the GNSS baselines. The weight matrix can be considered as a diagonal matrix, and the weight of each baseline k is defined as pk=1/σk2, where σ can be related to the length of baselines. Kuang [13, p. 183] suggests that the value can be between 0.1 and 2 ppm for deformation monitoring purposes. The E matrix is used to apply the inner constraint condition to the network. For more details about the structure of A, P and E matrices in networks of GNSS observations see Alizadeh-Khameneh et al. [1].

In the optimization procedure, the difference between the existing VC matrix of the network and the ideal criterion is sought to be minimum. In other words, the design problem tries to find an optimal configuration (A) and weights (P) such that the criterion matrix Cdc can be best approximated by Cd, i. e.,

where ∙2 stands for the L2-norm.

The contents of Cd in Eq. (22) are nonlinear functions of both the observation weights (pk) and coordinates of the network points (xi,yi,zi), which must be linearized by expanding to the Taylor series as below (Kuang [13], p. 267):

where Cd0 is the approximation of VC matrix that can be computed by Eq. (21) and considering an initial weight matrix; m and q are respectively, the number of points and baseline observations in a network. As already mentioned, we only deal with the second-order design problem in this study. Therefore, the terms related to the station coordinates are removed from Eq. (23) and only the last term of the equation, which shows the derivation of VC matrix with respect to weights (p) will remain.

An optimal design should also fulfill the reliability of the network and empower it in detecting and revealing possible gross errors. The reliability R is defined as the redundancy of observations as:

where I3m is an 3m×3m identity matrix.

The goal is to maximize the minimum reliability of the network, i. e.

where ∙∞ represents the L∞ norm.

To be able to embed the precision criterion matrix and the reliability and physical constraints in the optimization procedure, a bi-objective optimization model can be developed by considering the precision and reliability objective functions in Eqs (22) and (25) respectively (Kuang [13], p. 253).

Subjects to constraints

with

where u is a vector including the differences of the VC matrix and the criterion matrix in Eq. (22). The improvements of weights are stacked in the vector w. This vector is used to update the P matrix during the optimization procedure. The H matrix contains the derivatives of the Cd matrix with respect to the weights, pk. Using the vec∙ operator in the above equations stacks the columns of a matrix one below each other to build a vector. The Q1 and Q2 matrices provide and control the quality assurance of the network. H1 and u1 represent the precision control for the network; R2 is a matrix for accumulating the derivatives of the reliability matrix R with respect to the weights, R1 and ro are the vectors containing the diagonal elements of R and the possible minimum value for the reliability, respectively; A00 and b00 fulfill the physical conditions to avoid negative weights for estimated observation weights. The detailed contents of these matrices are described in Alizadeh-Khameneh et al. [1].

Equation (26) defines a bi-objective optimization model that is constrained to the precision and reliability and can be solved by applying a quadratic programming (see for instance, Milbert [16]). The first term of this equation represents the precision criterion and the second term shows the reliability criterion. A single-objective optimization model can be developed by using any of these terms individually. The advantage of different optimization models has been studied in, for example, Alizadeh-Khameneh et al. [1] and the bi-objective model is concluded to be the best model that can fulfills all the network quality criteria including the economic aspects.

6 Numerical studies

The concept of the finite element method can be used in geology and geodesy to analyze the strain parameters in deformation studies. However, adopting this method to design and optimize a geodetic network for deformation monitoring purpose is innovative and provides an effective solution to detect the deformation parameters in large-scale regions. In this study, we numerically implemented this solution in the Skåne area of Sweden, and designed an optimal network, with respect to precision and reliability.

6.1 Study area

The Skåne GPS network in Sweden, which previously was the case study for Pan, et al. [17] and Kiamehr & Sjöberg [12], is here chosen to assess the practicability of the method. The area is located in the southernmost part of Sweden and is expected to be one of the most seismically active fault zones in the country. Since 1989, the area was monitored epoch-wise by a GPS deformation network until 1998 [12]. The Skåne region, the location of seven established GPS stations (spaced approximately 80 km apart), and the approximate location of the fault zone – the so-called Tornquist zone – are depicted in Fig. 2. The data for this experiment was collected from two observation epochs of 1996 and 1998. The coordinates from these two epochs are available in ITRF96 as geocentric Cartesian coordinates, which are transformed to the Transverse Mercator map projection, used as a local coordinate system in this study. The height of the points in this local system is provided by the ellipsoidal heights.

Figure 2

The GPS monitoring network of Skåne in the South of Sweden. The right figure shows the location of established GPS stations and the fault zone in the area.

6.2 Results

The Skåne GPS deformation network is an example of regional deformation studies, where the monitoring procedure is carried out to detect the creeps between tectonic blocks along the faults. Usually, such monitoring network consists of a geodetic network, plus some isolated non-geodetic observables (for instance, measurements from strainmeters) that may not be geometrically connected with the geodetic network (Kuang [13], p. 181). However, geodetic measurements are adequate for deformation monitoring, and non-geodetic approaches can be supportive in providing some additional geological information about the deformable body. In the prior study of the Skåne deformation monitoring network by Pan, et al. [17] the displacement observations are studied with reference to a fixed point in Onsala, Sweden, located about 180 km away from the area, while in this study, we treated the network as a relative network, where all the GPS observation points are located in the area. Therefore, the whole area is subject to deformation. This assumption enables us to divide the whole area into small blocks and estimate the deformation components in each element based on the explained methodology in Eqs (8) to (16).

Figure 3

The Skåne area, which is discretized to ten subareas with the help of the 3D Delaunay triangulation technique.

One of the imperative and important steps in optimizing a geodetic network is to define an appropriate precision criterion. For this purpose, we first need to gain some insight into the magnitude of deformation parameters in the deformable body. Therefore, we started by splitting the body to be able to implement the FESA in each element and consequently obtain the required information. As a first step, the Skåne network is split into a finite number of optimal non-overlapping elements. With the help of the 3D Delaunay triangulation, the area is split into 10 elements, where the network points are considered as the vertices of the elements. The Delaunay triangulation method expresses which four points compose an element. Figure 3 shows the region divided into a number of sub-regions. For a better illustration of the discretized area and avoiding a cluttered image, only one of the tetrahedral elements is shown with a solid color in the figure.

Table 1 provides information on the output of the Delaunay triangulation. The number of finite elements and the network points that are involved in constructing them are shown in this table.

Table 1

The number of tetrahedral elements by the Delaunay triangulation.

Element No.	Point Connection	Point Names (ref. Figure 3)
I	2-6-3-4	1: Stavershult 2: Limhamn 3: Kivik 4: Dalby 5: Veberöd 6: Solvik 7: Hässleholm
II	7-3-5-4
III	6-2-1-5
IV	5-2-1-4
V	2-3-7-4
VI	3-6-5-4
VII	6-2-5-4
VIII	7-5-1-4
IX	3-7-1-4
X	3-5-1-7

Table 2

Estimated deformation parameters and computed dilation values for the 10 split elements of the Skåne area.

Element No.		I	II	III	IV	V	VI	VII	VIII	IX	X
Strain Parameters	ex	2e-7	1e-7	-2e-6	-9e-7	2e-7	2e-7	-8e-7	4e-7	-6e-8	2e-8
	ey	-3e-9	2e-7	-4e-8	-2e-7	4e-8	-4e-7	-2e-7	5e-8	1e-7	2e-8
	ez	1e-3	1e-3	6e-3	1e-3	6e-4	7e-4	2e-4	9e-4	5e-4	3e-4
	exy	-7e-8	1e-7	-4e-7	-5e-7	9e-8	-3e-7	-6e-7	2e-7	8e-8	6e-9
	eyz	-3e-5	7e-5	3e-5	9e-5	-3e-5	3e-5	1e-4	5e-5	-1e-5	9e-6
	exz	-5e-5	9e-5	2e-4	1e-4	-7e-5	2e-5	1e-4	6e-5	-3e-5	7e-7
Differential Rotations	ωx	-3e-5	7e-5	4e-5	9e-5	-3e-5	3e-5	1e-4	5e-5	-1e-5	9e-6
	ωy	-5e-5	9e-5	2e-4	1e-4	-7e-5	2e-5	1e-4	5e-5	-3e-5	2e-6
	ωz	-2e-7	7e-8	-3e-7	1e-7	-5e-8	-4e-7	1e-7	-1e-7	1e-7	4e-9
Dilation	Δ	1e-7	1e-7	-9e-7	-6e-7	1e-7	-1e-7	-5e-7	2e-7	4e-8	2e-8

Table 3

The standard deviation of deformation parameters for the 10 split elements.

Element No.		I	II	III	IV	V	VI	VII	VIII	IX	X
Standard deviation of deformation components	σεx	1e-13	7e-14	1e-12	4e-13	9e-14	8e-14	6e-13	2e-13	9e-14	1e-13
	σεy	2e-13	1e-13	4e-13	2e-13	7e-14	2e-13	2e-13	5e-14	8e-14	6e-14
	σεz	7e-11	9e-11	3e-10	1e-10	6e-11	5e-11	2e-10	5e-11	4e-11	5e-11
	σεxy	1e-13	7e-14	6e-13	2e-13	6e-14	1e-13	3e-13	1e-13	6e-14	6e-14
	σεyz	4e-11	5e-11	1e-10	6e-11	3e-11	2e-11	8e-11	3e-11	2e-11	2e-11
	σεxz	4e-11	5e-11	1e-10	6e-11	3e-11	2e-11	8e-11	3e-11	2e-11	2e-11
	σωx	4e-11	5e-11	1e-10	6e-11	3e-11	2e-11	8e-11	3e-11	2e-11	2e-11
	σωy	4e-11	5e-11	1e-10	6e-11	3e-11	2e-11	8e-11	3e-11	2e-11	2e-11
	σωz	1e-13	7e-14	6e-13	2e-13	6e-14	1e-13	3e-13	1e-13	6e-14	6e-14

The estimated deformation parameters and their corresponding standard deviations are presented in Table 2 and 3. The strain parameters in these tables do not have units and can be expressed as a relative unit. The unit of differential rotations is radians. The aim of this paper is not to analyze the magnitude of the deformation in this area, and therefore no detailed interpretations for the obtained results are provided here. For more information on the geological investigation of the area see Pan, et al. [17]. However, a deformation primitive – dilation – can be derived from the strain matrix to analyze the average expansion/contraction (+/-) in each element and obtain an overall information on the deformation of the area.

The dilation parameter, which is provided in the last row of Table 2, shows that the minimum deformation occurs in the elements IX and X. Recalling Table 1 and Fig. 2, it can be seen that these elements belong to the upper part of the Tornquist zone, and this is in agreement with previous deformation analysis of this area by Pan, et al. [17]. They concluded that the southern part of the Tornquist zone is more exposed to deformation than the northern part. According to Table 2, the major deformations belong to the elements III and IV, located in the southern part.

To design a deformation network, it is required to acquire some a priori geological information on the area of interest besides the other geodetic information. Notwithstanding lack of a priori geological information in this study, the required information is obtained from the two-epoch geodetic measurements of the network – presented in Table 2 and 3. According to Table 3, an average optimization boundary for the precision of deformation parameters could have been defined similar to, for instance, Kuang [13], p. 294, where a standard deviation of 1 ppm is chosen for all deformation parameters. However, the innovation in this study is to use different precision boundaries for each element, thus we assumed the values in Table 3 as the precision boundary values for the optimization procedure and used them in forming the VC matrix of deformation parameters Ce in Eq. (17). This enables us to define a criterion matrix Cdc, where different displacement precision is sought on each network point according to a priori deformation information of its adjacent elements.

Moreover, it is also possible to obtain the VC matrix of deformation parameters matrix, i. e. Ce, based on the standard deviations of the stress parameters, as the stress is an exerted force to a body that can cause a change and deformation in its shape. This can be achieved by using Hooke’s formula and applying Young’s modulus and Poisson’s ratio [15]. Young’s modulus is a measure of the stiffness of a solid body, thus to obtain its relevant value and implement it in the network design problem, it is needed to have some a priori information on the area.

The precision of the displacements is defined based on the precision of the deformation parameters in each of the ten-split elements to be the precision criterion of the network. The reliability threshold is also set to 0.65, meaning that the reliability of each baseline should not be less than this value The optimization procedure maximizes the minimum reliabilities and increases their values (Eq. (25)). The reliability of an observation varies in the range of 0 to 1 (minimum to maximum), and the mentioned threshold is an empirical choice for this study. It is also possible to decrease the reliability of other observations, but the procedure keeps their value above the threshold. As the second-order design stage is performed here, the optimum baseline weights are obtained after the optimization. Since every baseline in GNSS measurements has three components, its corresponding weight matrix consists of three weights for each component. Therefore, in the optimization procedure, it is of importance to remove all components of the weight matrix for unnecessary observations.

Table 4

Comparison between different optimization models, i. e. Single-Objective Optimization Model (SOOM) and Bi-Objective Optimization Model (BOOM).

Objective Function	fprecision	freliability	No. of Removed baselines out of total 21 baselines
SOOM of Precision	8.34	1.16	9
SOOM of Reliability	8.53	1.95	1
BOOM of Precision and Reliability	8.42	1.41	4

Three different optimization models are applied to the network: two single-objective optimization models of precision and reliability, and one bi-objective model of both precision and reliability. The bi-objective model is more efficient as it overcomes the inconsistency between the precision and reliability criteria and it simultaneously fulfills both quality criteria of the network. The values of fprecision (objective function for precision) in Table 4 represent the norm of differences between the network existing precision and the criterion. As expected, the single model of precision achieves the best fit of these two matrices. Similarly, the norm of the reliability vectors is shown in the table under freliability (objective function for reliability). It can be seen that the higher reliability value belongs to the single model of the reliability. However, the bi-objective model keeps a balance and fulfills both criteria. On the other hand, the number of eliminated baselines is of importance from an economic point of view. Figure 4 illustrates the optimized observation plans of the Skåne network with the number of baselines eliminated after applying the three different optimization models. The reliability model (plot b) is the strictest one and retains almost all the baselines – excluding only one – to ensure the capability of the network in detecting errors. It seems that the bi-objective model (plot c) keeps the balance here as well and removes unnecessary baselines more than the reliability model but less than the precision model.

Figure 4

The observation plan of the Skåne network after applying the three different optimization models. The red dashed lines represent the eliminated baselines from each observation plan, and the blue solid lines are the necessary baselines to be observed. The used optimization models in the plots a to c are: SOOM of precision, SOOM of reliability, and BOOM of precision and reliability.

According to the defined precision criterion in Eq. (17), the displacement of the network points can be measured with the precision of minimum about 3 mm in this network designed by the bi-objective optimization model. In Fig. 5, three graphs are provided to show the displacement precision in the network before and after the optimization. The criterion does not have similar values for all the points as it has been defined by considering the precision of the deformation parameters in the different 3D elements. In other words, the precision criterion of the displacements at each point in this study is under influence of the precision of the deformation parameters of the adjacent tetrahedral elements.

Figure 5

The precision of network points before and after performing the bi-objective optimization model – the blue squares and green circles, respectively. The red stars represent the precision criterion. Network point ID numbers refer to Table 1.

It can be seen in Fig. 5 that the lowest precision belongs to the point number 4. Referring to Table 1, where the dissected elements are shown, we can see that this point is involved in the construction of most of the elements. Thus performing the precise measurement from this point is essential for fulfilling the precision of deformation parameters in different elements. Now, if we go back to Fig. 4, we can see that when the objective function of the optimization model is only the precision (plot a), the optimization procedure retains the all observations from the point number 4, as a tie point for all the dissected finite elements. Moreover, the precision criterion for the elements III and IV, which have larger deformations (according to Table 2), are determined to be higher.

7 Conclusion

An optimal design of a geodetic network helps the surveyors to establish a network, which can fulfill the quality requirements of the network. Furthermore, the optimization procedure suggests an observation plan that has a lower cost. Dealing with deformation monitoring networks, an optimally designed network may save a considerable amount of cost and energy in the project as the observations in these cases must be repeated in continuous or epoch-wise measurements.

Using geodetic measurements, it is possible to measure the displacement of discrete network points in different time intervals, which are occurred due to deformation of a deformable body. To estimate the deformation parameters (strain and rotation), a proper deformation model must be fitted to the displacement field of the object. A displacement field is a difference between the two body states. In case of large-scale areas (regional or global networks), it is expected to observe different crustal movement behavior within an object. Thus, it is reasonable to discretize the region of study into sub-regions (finite elements). One of the well-known approaches for this purpose is the Delaunay triangulation, which can split the deformable body into a number of triangles or tetrahedrons. We used this technique in this paper to define a precision criterion matrix in optimizing the Skåne GPS monitoring network. The area of the network was split into 10 finite elements with the help of the Delaunay triangulation.

Using the precision of deformation parameters directly in the optimization procedure is possible, when the whole area is assumed as one deformable body. However, this solution can be very complicated, when the area is split into a number of elements. Therefore, a methodology is developed to define a criterion for the precision of displacements of the network points based on the required precisions of deformation parameters in each of tetrahedral finite elements.

The defined criterion here enabled the network to be sensitive for minimum 3 mm displacement at each point. A bi-objective optimization model is then developed based on the precision and reliability criteria to optimize the Skåne GPS network. The initial number of 21 GPS baselines was decreased to 17 after performing the optimization procedure. It should be mentioned that there are only seven network points distributed in the whole area. It is obvious that considering more GPS stations within the area yields more baseline observations and potentially more savings after performing the optimization procedure.