Hybrid algorithm for identification of stable reference points in a monitoring network

1 Introduction

Many reasons cause the earth’s surface deformations. They can be natural or indicated by human activity. Whatever the reason, the changes can affect technical infrastructure, cultural heritage, or in some cases, pose a danger to human life. These are the reasons for monitoring the deformations to evaluate safety level and predict potential hazards. Water dams [1], [2], [3], bridges [4], 5], industrial structures [6], tall buildings [7], and cultural heritage sites [8], 9] are the most common objects of deformation monitoring. Deformation monitoring can also be applied to better understand natural processes like Earth’s crust movements [10], landslides [11], or glacier changes [12], especially when connected with the safety of engineering structures [13], 14].

A wide variety of measuring tools and measuring methods can be applied for deformation monitoring. Some of them are commonly classified as geotechnical (structural) and include sensors such as accelerometers, extensometers, inclinometers, strainmeters, and others [15], 16]. These tools are able to register physical quantities as a result of deformations. They usually provide information about local changes and are often applied in automatically working control systems monitoring the object’s health in real time.

However, there are situations where the absolute displacements are necessary. In such cases, geodetic methods can be applied. Geodetic measurements allow for the registration of the geometrical state of the measuring network, which consists of control points on the object and reference points that are expected to be stable. Measurement of the network is repeated periodically. The measurement of the network completed at a specific time is referred to as the measurement epoch. The object deformations can be obtained by comparison of two or more epochs. The object points displacements can be considered absolute only unless a part of the reference points remains stable between two analysed epochs.

Although the network design assumes that the location of the reference points is free of deformations, this assumption is often verified negatively in real situations. Such a situation is more likely when a long period is considered. Therefore, verifying reference points’ stability appears as a crucial element preceding the calculation of displacements. Many methods of stability analysis have been developed [17]. Identification of the stable reference points can be conducted by determining the displacements, e.g., [18] or separately, as a procedure preceding the calculation of displacements. Some identification algorithms use observation values of analysed epochs [18], [19], [20], and others are based on coordinates, which usually come from an independent preliminary adjustment of individual epochs. In the last case, the precision of the coordinates is usually available, and the coordinate transformation procedure is used for stability analysis [21], [22], [23].

Regardless of the applied identification algorithm, the problem is not simple, and the identification result can be equivocal. This may occur when the network covers a wide, particularly geologically differentiated, area [24], 25]. Additional source of difficulties may be varied in the direction character of displacement causes. As a result, more than one group of mutually stable but unstable concerning other groups, reference points can be identified. The commonly applied identification algorithms usually do not deal properly with such cases.

This article presents a new algorithm that is based on coordinate transformation and is able to detect various groups of mutually stable points. It is a hybrid algorithm, and it uses optimization methods to search for the solution. The proposed algorithm is described in Section 2.2.3. The key element that enables the application of optimization algorithms is the objective function. Two variants of the objective function were proposed in this article. They are described in Section 2.1 and they are both based on the weight functions originally applied for the robust estimation.

The proposed algorithm is applicable for vertical (1D), horizontal (2D), or spatial (3D) networks, but due to the simplicity of the vertical displacements (1D) measurements, it is reasonable to consider its application for 2D or 3D networks.

The detailed description of the proposed algorithm is presented in Section 2. Section 3 shows the results of the numerical test for the simulated 2D object, which was conducted to verify the algorithm’s effectiveness. The discussion of the results and the conclusions are presented in Section 4.

2 Algorithm of stability analysis

Measurement results of each epoch are the subject of preliminary adjustment. This process verifies the consistency of internal measurement data and delivers the coordinates of the object points along with their accuracy assessment. Local network reference (e.g., free condition) is usually applied to the preliminary adjustment.

The results of the preliminary adjustment can be used for the deformation analysis. In this approach, the coordinate differences of the control points represent their displacements. This is true only when the coordinates of two analysed epochs are referred to the same stable reference system. It can be achieved by the coordinate transformation using stable reference points as the transformation tie points. Coordinate transformation can be defined by the following equation:

Where:

a _i – position vector for the point in the starting system of coordinates,

b _i – position vector for the point in the final system,

t – translation vector.

R is a rotation matrix that is applied only to 2D and 3D transformations.

Assuming that vector a _i in equation (1) corresponds with the coordinates of point i for the first epoch, and vector c _i with coordinates for the same point for the second epoch, residuum vector r _i (2) can be determined for each point:

The residual vector r _i is the indicator of the control points’ stability. If there exists a group of control points that maintain mutual consistency in both considered epochs, it is possible to find such transformation parameters that residuals r _i for them will not excide values resulting from measurement precision. The final task of the reference system identification is to find such transformation parameters represented by R and t that the residual vectors for the most numerous group of potentially stable reference points will not exceed the assumed threshold:

where ||r _i || is the L2 (Euclidian) norm of residual vector. As the point coordinates for both epochs come from preliminary adjustment, the positional precision can be used to formulate the threshold value:

The factor m denotes the relation between the precision characteristics resulting from preliminary adjustments of both epochs and the threshold connected with the expectation of the points’ stability. In the model case, when the group of stable points is not affected by any displacements, the value of m is often taken as 2. In practical cases, especially when the analysed epochs are divided over a long period of time or when the terrain is not stable over a large area, the criterion of stability, except for the network measurement precision, must also take into account other influences as the control points’ micro-movements. The value of m has to be enlarged then. Such a case is described by García-Asenjo et al. in [14].

Knowing the final transformation parameters by which residua on reference points are within an acceptable range, we can easily calculate the displacements of all control points as transformation residuals r i = r X i , r Y i T (2). The spatial movement length dp _i can be determined from equation (6).

There are several ways to obtain the desirable transformation parameters. The first one is a commonly known congruency test [26], [27], [28]. Other methods that can be listed are numerical control of object adherence [29], Iterative Weighted Similarity Transformation (IWST) [5], 30], 31], or a robust algorithm proposed by Hu et al. [32], using Mahalanobis distance.

Another category that the method proposed in this article belongs to includes methods using optimization algorithms. They are based on the assumption that transformation parameters create a certain space of solutions (search space). The dimension of such space depends on the transformation type, and it can be one for 1D transformation, three for 2D transformation, or six for 3D transformation. Each point of this space corresponds to a set of transformation parameters and, consequently, to the set of residuals on all considered tie points. In the work [23], two optimization algorithms were proposed – the Simulated Annealing and Hooke–Jeeves algorithm. Both algorithms used a new objective function that considered values of transformation residuals and the number of congruent points.

Another issue that must be taken into account is the potential existence of various mutually stable groups of control points. It can appear when the network covers a large, differentiated in geological structure area. Such cases are the subject of [24] and are especially challenging for identification procedures. Selecting the dominating, most numerous group can be incorrect in specific situations. In the author’s opinion, only an expert can make the right decision, taking into account aspects beyond the geometry and precision of the network.

The desired identification algorithm should, therefore, detect all mutually stable subgroups with the transformation residuals for each subgroup. Such an approach was applied in [23] and is continued in this study.

2.1 Objective function

In the case of coordinate transformation, an objective function delivers a measure of the fitting of congruent points, which is represented by the transformation residuals. The definition of the objective function result is a complex question. The measure of ‘fitting’ must take not only the set of values of transformation residuals but also consider that residuals for a group of unstable points (the largest) should not affect the measure of the fitting of the congruent points.

In [23], the author proposed an objective function based on the L1 norm of coordinate transformation residuals. If a group of congruent points is detected, some points with the maximal residuals are excluded from the L1 norm. The function appeared to be efficient in detecting the groups of congruent points, but it has two flaws. Selection of the maximal residuals needs ordering the points by the length of the residual vector, which increases numerical cost – especially when the function is used many thousands of times, as it is in the case of optimization methods. Furthermore, the rejection of the large residual vector from the objective function causes a sudden change in the objective function value, which hinders the operation of some optimization procedures.

The mentioned disadvantages caused the new research to be undertaken, aimed at finding a new, improved, objective function formulation. Due to the problem specificity, functions applied in the robust estimation were taken into consideration. The robust estimation is a method of network adjustment with the assumption of the existence of outlying observations. The estimation process is conducted iteratively, and the basic element of the adjustment procedure is a weight function that modifies the observations’ weights depending on their residuals. If the observation gets a large residuum, it gets, in consequence, decreased weight, and it does not influence the estimation results in the next iteration. There are many robust estimation methods, such as Huber’s [33], Hampel’s [34], Danish [35], Kadaj’s [36], and many others. They differ in the weight function definition, the effectiveness of the parameters estimation, and outliers elimination.

Applying the robust approach for deformation analysis is not new. Except for observation adjustment [36], [37], [38], they were used for establishing the deformation model [3], 39], 40].

Two robust functions have been applied in this research. The first one is the classical Huber weight function. It can be expressed by formula (7),

(7) h x = 1   if ‖ x ‖ < f f ‖ x ‖   otherwise

where x is the residual vector, and f originally stands for an observation’s standard deviation, and, in this study, is connected with preliminary stability assessment. The Huber weight function is commonly used in many applications. It is quite simple, but it is composed of two elementary functions that require the check of the argument value to choose the corresponding elementary function.

Equation (8) describes the Kadaj weight function [36], where k _x originally denotes the standard deviation of x, and in this study, it is connected with a preliminary assessment of stability. Although the form of the function is more complex (it contains the power function), it is uniform in the whole range of argument values.

(8) k x = e − x 2 2 k x 2

The objective function is expected to show the places of the search space that correspond with the objective, which, in most solutions, is associated with the minimal value of the function.

The selected robust weight functions originally have a maximum of x = 0. To obtain a minimum for the small residuals, the original functions must be modified, as shown in equations (9) and (10). To compose the objective function for which a minimum is associated with the objective, we must formulate auxiliary functions H(x) and K(x) by reversal of the original robust weight functions.

(9) H x = 1 − h x = 0   if ‖ x ‖ < f 1 − f ‖ x ‖   otherwise

(10) K x = 1 − k x = 1 − e − x 2 2 k x 2

The distribution of both auxiliary functions is shown in Figure 1.

Figure 1:

Distribution of auxiliary functions for f = 1 and k _x = 2.

Assuming r _i is the vector of transformation residuals the corresponding objective function for the Huber weight function can be formulated as:

(11) F H = ∑ i = 1 n H ‖ r i ‖

The detailed analysis showed that the Kadaj function (10) reaches its asymptote very fast (for x/k = 9, K x = 1 with numerical precision of 10⁻¹⁷). It is not a problem in the case of robust network adjustment, but when applying for constant points identification, it does not allow for distinguishing between the solutions that are far from the optimum. This fact can cause a problem for optimization algorithms that use the comparison of objective function values to approach the solution. To overcome this issue, an additional element of the L2 norm was introduced into the objective function based on the Kadaj weight function (12). As a result, the objective function F _K gets different values even for large r, for which K ‖ r i ‖ = 1 .

(12) F K = ∑ i = 1 n K ‖ r i ‖ + c ‖ r i ‖ 2

The value of the c parameter is selected empirically. It should be as small as possible to distinguish various solutions far from the optimum and not to disrupt the main objective function distribution.

2.2 Optimization methods

Optimization methods are applied in many areas of science and technology. They allow for finding the best solution when the problem has its mathematical model defined by certain number of parameters. The parameters create the space of solutions. Each point of this space is described by the objective function value, indicating the level of reaching the objective. A wide variety of algorithms have been invented to look for the solution. Their effectiveness mainly depends on the number of parameters and objective function distribution. The rapid development in computational technology induced the development of optimization, resulting in new algorithms, often inspired by natural processes like evolutionary rules [41], 42] or living organisms’ behavior [43].

Considering the number of parameters, the coordinate transformation problem is not very complex. The maximum number of parameters in the case of 3D transformation does not exceed 7. It allows for the application of some algorithms that lose their effectiveness when the number of parameters is increased. The objective function value distribution is difficult to predict. It can have one or more minima. This suggests a choice of algorithms of universal character.

Two algorithms were applied in [23] for the solution of constant points identification. The Hooke–Jeeves and Simulated Annealing algorithm. Although simulated annealing appeared to be more effective in the analysed problem solution, the Hooke–Jeeves algorithm also appeared to have some advantages. Both algorithms are described in the subsequent Sections 2.2.1. and 2.2.2. A combination of both algorithms using the main advantages of each of them is proposed in Section 2.2.3.

2.2.1 Hooke–Jeeves algorithm

The Hooke–Jeeves (HJ) local search algorithm proposed in 1961 [44] is an efficient way to find the optimization problem solution. The HJ is referred to as the pattern search algorithm, and it is fully deterministic. It is simple and relatively fast-converging. Because no objective function gradient is necessary for the algorithm to work, it can be convenient for cases when the objective function has no analytical form and is obtained based on empirical data.

The following parameters must be defined to apply the HJ method to a particular problem:

d – orthogonal basis of n independent linearly orthogonal vectors,

τ₀ – initial length of the search step, dependent on the search space as well as the objective function distribution,

γ – search step decreasing factor in subsequent iterations,

τ_end – minimal step length constituting the criterion for ending the iterative process,

x₀ – the starting point of the procedure.

A more detailed description of the algorithm can be found in [23] or [45]. The HJ method is commonly used to solve optimization problems in various fields of engineering. Identification of the water dam parameters [46] or aircraft parameter estimation [47] can be given as examples. Applications for the problems of geodetic measurements are not numerous, but 3D transformation problems [48], identifying reference points in displacement monitoring networks [23], or second order design of geodetic networks [45] can be quoted as examples.

2.2.2 Simulated annealing

The algorithm was first proposed in the work of Metropolis et al. [49], but only the development of computational technology made it possible to apply the algorithm to solve real problems. The work of Kirkpatrick et al. [50] was also an important milestone. The simulated annealing (SA) belongs to the metaheuristic algorithms. Its main advantage is the ability to find the global minimum when local minima exist. The SA algorithm is also useful when the solution has to meet complex criteria.

Its features caused the SA method to be applied to solve a wide variety of problems from various areas of technology and science. Solving a system of multiple linear equations [41], structural health monitoring [51], a bearing fault diagnosis [52], or UAV path planning [53] are only selected examples from a large number of others. The algorithm is also used for the solution of various problems in the area of geodetic measurements, like geodetic measurements results procesing [54], design of measurement networks [43], 45], 55], 56], deformation analysis [23], or coordinate transformation [57]. A thorough description of SA is provided in [58].

The basic idea of the SA algorithm is iterative searching through the search space, looking for the best solution reflected by the optimal (usually minimal) value of the objective function. The new, potentially better solution is generated randomly using normal distribution centered in the current solution. The search scope corresponding to the temperature decreases with subsequent iterations. Assuming x _i is a current solution, Δx _i is a randomly generated incremental vector to the potential new solution the relation (13) describes the way of choosing a new solution.

(13) x i + 1 = x i + Δ x i   if f x i + Δ x i < f x i x i + Δ x i   if f x i + Δ x i ≥ f x i with probability p x i   otherwise

The new solution is generally the solution corresponding to the lower value of the objective function (case 1). The algorithm allows, however, the possibility of choosing a worse solution (case 2). It can be applied in an extremely complex distribution of the objective function with many local minima when the global minimum is the only aim of the process. The probability p of the choice of the worse solution is determined with the Boltzmann distribution. In simplified cases, a small constant value e.g., p = 0.001 or p = 0 can be assumed.

Application of the SA method to the stability analysis requires defining a few key elements.

2.2.2.1 Definition of the objective function

The objective function evaluates the current solution, allowing the approach to the solution. In this work, the objective function is described in Section 2.1.

2.2.2.2 Cooling scheme

The cooling scheme defines the initial temperature and the formula for determining the decrease in temperature in subsequent iterations. The separate factor t is assumed to be a temperature in this application. t₀ = 1 was taken as the initial value of t. The change in t values in subsequent iterations is described with formula (14),

(14) t i = t i − 1 β = t 0 β i

where i denotes the iteration number, and β is the so-called cooling factor. Its value is selected depending on the nature of the optimization problem.

2.2.2.3 Molecule movement model

Two elements must be defined to define the molecule movement model. The first is the way of generating a new solution. As it is commonly accepted, a new solution is generated in the way described by formula (15):

(15) x i + 1 = x i + Δ x i

In this study, the particular transformation parameters constitute the vector x. The incremental vector Δx _i is generated randomly, using normal distribution N(0, σ _i ), where σ _i is the vector of standard deviations of particular transformation parameters corresponding to x in i-th iteration. Formula (16) describes the value of the standard deviation in iteration i:

(16) σ i = σ 0 t i = σ 0 t 0 β i

The initial values of standard deviation σ₀ are selected depending on the nature of the problem, and they reflect the assumed range of the corresponding parameter variability.

The acceptance criterion for a new solution is the second element of the molecule movement model. In this study, formula (13) was applied with p = 0, which means there is no possibility of accepting a worse solution.

2.2.2.4 Criterion for ending the iterative process

The condition for ending the iterative process can be formulated based on:

1. the final temperature,

2. the number of iterations,

3. the acceptable value of the objective function

In this study, the temperature threshold t_end is applied (17).

(17) t i < t end

The SA algorithm is able to find a global solution when local solutions exist. However, in this study, the local minima corresponding to the subgroups are also an object of interest. The ability of the SA algorithm to find local minima except for global minimum can be controlled by the cooling speed. Slowing down the cooling process (colling factor β closer to 1) makes the procedure more likely to find a global minimum. In case of faster cooling, the procedure will more likely return the solutions (local minima) corresponding to the smaller subgroups of mutually stable control points.

In conclusion, there is a challenging task to select the optimal value of the cooling factor to detect the real subgroups of the stable control points and minimize the junk solutions that do not correspond with the acceptable subgroups.

2.2.3 Hybrid algorithm

The optimization methods vary in terms of algorithms used, which are dedicated to the specific types of problems. It also applies to the two above-described algorithms. The SA needs a large amount of calculations. It is able to find the optimal (global) solution with a complex distribution of the objective function and when the local minima are present. On the other hand, it is not very effective when the objective function distribution is regular with one evident minimum. The second algorithm (HJ) is considered fast, but it is likely to find a local minimum as a final solution.

In the case of the stable points identification with the objective function described in Section 2.1 multiple local minima are possible. They can correspond to the separate subgroups of points that are mutually stable within a group but are not stable in relation to the points of another group. Also the single unstable points that do not belong to any group can create a local minima in the search space of transformation parameters. Such situations can appear in real cases [24], 25].

The existence of multiple subgroups can be problematic for some traditional algorithms based on the assumption of existing one group of stable points that contains the majority of the reference network. The example concerning the commonly known congruency test was presented in [23]. The algorithms proposed in the works of Neitzel [20], 59] or Lehman and Lösler [60] are more efficient. They can identify even very small subgroups of congruent points, but they use a combinatorial approach, which leads to an exponential increase in computation cost with the number of points. In consequence, the combinatorial methods can be applied only for not very large networks. In the case of optimization algorithms, the computation cost also depends on the network size, but this dependence is only linear.

As only one solution (group of congruent reference points) can be found in the single activation of the hybrid procedure, it must be activated a number of times to properly identify all the potential, internally congruent groups of reference points. This number depends on the existence of smaller subgroups that are identified with lower probability. Based on the series of experiments, a number of trials of 500–1,000 can be recommended.

The research [23] shows that although the minima corresponding with the particular groups are evident, the distribution of the objective function in the minimum vicinity is more complex and may contain local minima. As a result, the HJ algorithm delivers quite a large number of solutions corresponding to the small subsets of the larger groups of mutually stable points or with unstable points. The SA algorithm performs better, but it suffers from a far from optimal procedure of selecting the next step at the initial stage of the search process, when the way to the closest minimum is obvious. Another problem connected with the SA algorithm is its tendency to find global or dominating minima with a lower objective function value. In the case of this study, where all possible subgroups of stable points should be found, it is not an advantage.

That was the reason for proposing the new hybrid algorithm, which is a combination of two previously described algorithms. The idea of using hybrid algorithms is not new. Twardochleb et al. [61] proposed a combination of Monte Carlo and HJ methods to solve problems of an economy. Mustafa [62] applied a combination of two metaheuristic algorithms – Particle Swarm Optimization and SA for the maximum partition problem. Rios-Coelho A.C. et al. [63] applied a combination of Particle Collision Algorithm with HJ for the chemical equilibrium nonlinear systems. Those and many other examples show that a combination of appropriately selected elementary algorithms can increase the efficiency of the problem solution.

In the case of this study, the combination of HJ and SA algorithms was used to establish a new hybrid algorithm. The HJ is the main method. It is used to find the first (closest) minimum of the objective function, which is then a starting point for the SA algorithm. The SA algorithm plays a supporting role and was modified for application within a hybrid algorithm. It starts at the minimum point found by HJ, and its function is to test whether the minimum is not a subminimum. The modification concerns the termination criteria. The algorithm finishes its work in case a temperature reaches the assumed threshold (minimum confirmed) or when a point with a lower value of the objective function than the current minimum is found. In the second case, the point found becomes the starting point from which the subsequent run of HJ is started. As the SA algorithm is not intended to find the final solution, its cooling rate can be faster, which shortens the whole procedure. The flowchart of the modified version of the SA algorithm is shown in Figure 2. The hybrid procedure finishes when SA is not able to find the new minimum and thus confirms the minimum found by HJ. The detailed flowchart of the proposed hybrid algorithm is shown in Figure 3.

Figure 2:

Flowchart of the modified SA algorithm.

Figure 3:

Flowchart of the hybrid algorithm.

The proposed hybrid algorithm has two main advantages. The first one is the shorter time of operation, which arises from the fact that the faster HJ method is applied to seek the solution, and a faster cooling rate is applied for the SA algorithm. Another advantage resulting from the features of the HJ algorithm is a higher probability of detecting the non-dominating subgroups of stable points that are also an object of interest.

As multiple minima can appear in the objective function distribution, it is necessary to set a strategy for how to distinguish the local minima corresponding to the searched groups of constant points from the others that are connected with smaller groups or single unstable points. Two parameters are of key importance to adjust the hybrid algorithm to the particular task.

The first one is the objective function parameter connected with the expected constant points congruency (k for F _K or f for F _H ). The first estimation of this parameter is connected with a preliminary assessment of stability. The experiments show that k = mσ _XY or correspondingly f = 0.5mσ _XY for m = 2 are the recomendable initial values. As it was mentioned in Section 2, in many cases, it is impossible to assess the stability level using the preliminary adjustment, and the objective function parameter must be adjusted empirically by gradually increasing the m parameter. There is no strictly formulated limit for m, but one must remember that the bigger the m value, the worse the final precision of the displacements.

Another parameter that strongly influences the hybrid algorithm’s performance is the cooling factor of the SA algorithm, which determines the ability to skip the local minimum found by the HJ algorithm. The effectiveness of the SA procedure depends, however, on the number, location, and stability level of the constant points. Due to the unknown number and location of the constant points, the cooling factor must be adjusted experimentally. The reasonable range of cooling factor values can be taken as 0.999 to 0.9999. In justified cases, the cooling can be even slower.

3 Numerical example

The proposed hybrid algorithm was tested on a numerical example. The simulated two-epoch 2D test network was used for the test. It contains 20 points. The points’ placement is shown in Figure 4. Table 1 contains point coordinates generated for both epochs. The network contains two subgroups of internally stable points. Points creating a particular subgroup keep mutual stability within the group, but are not stable in relation to the other points. 11 points are unstable in relation to both subgroups of stable points as well as to themselves mutually. The movements assumed for the groups are listed in Table 2. The coordinates were noised by true errors with the normal distribution and standard deviation of σ _XY = 5 mm, which is an equivalent of disturbance in both epochs with errors of 3.5 mm standard deviation. 2D transformation parameters from epoch 1 to epoch 2 were assumed as the solution vector elements. With the assumption of the constant scale, the solution vector x consists of three elements: x₁ = dX, x₂ = dY, x₃ = dα. Where dX and dY denote the incremental shift of the object’s barycenter in X and Y direction and dα denotes the rotation change between epoch 1 and 2.

Figure 4:

Test network points placement.

Figure 5 shows the distribution of the objective function for F _K objective function (a) and F _H objective function (b). Rotation angle α = −0.0040 g was assumed for the pictures. For better legibility, the values of both functions are inverted (the top of the hill in the picture corresponds to the minimum of the function).

Figure 5:

Resulting objective function. (a) – F _K , (b) – F _H

The experiment was aimed at checking if the proposed algorithm is able to effectively identify both subgroups of constant points. It consisted of performing n = 1,000 trials for three algorithms: HJ, SA, and a new hybrid algorithm. The results are shown in Tables 3, 4, and 5. Each algorithm was tested using both analysed objective functions – F _K and F _H . For the F _K  objective function (12), c = 1.0/m² was set. This choice ensures that the L2 norm component does not significantly affect the value of F _K when the components of r _i are expressed in meters. The frequency of finding, as well as the corresponding objective function value, is listed for each solution. As the same precision was assumed for all points, the common threshold value ɛ _i = 0.008 m was applied for all points and all tested algorithms.

3.3 Hybrid algorithm

The parameters for the HJ algorithm were the same as those for HJ as an individual algorithm analysis.

Similarly to the SA as a stand-alone algorithm, the same search space extent parameters were applied. σ 0 = σ x 0 , σ y 0 , σ α 0 = [0.15 m, 0.15 m, 0.015 g] was taken for both objective functions.

The value of the cooling factor was selected in an experimental way. For F _K , the cooling factor β = 0.9995 must be set for the proper identification of both subgroups. For F _H , β = 0.999 appeared to be enough.

The identification results of the Hybrid algorithm are summarized in Table 5.

As it can be realized, both variants of the objective functions based on robust estimation weight functions can efficiently show the places in the search space that correspond to the groups of stable reference points. Application of the F _K objective function means a higher computational cost, but the function distribution is more regular. The F _H objective function produces more local minima, which can be seen in Figure 5b, and leads to numerous solutions when using the HJ algorithm.

The test confirmed that the HJ algorithm is prone to irregular objective function distribution, which is especially visible in the case of the F _H objective function. Although many of the erroneous solutions can be filtered using quite simple criteria, it must be concluded that the HJ algorithm is not advisable as a stand-alone method.

The SA algorithm appears to be able to perform much better, but this was at the expense of longer computational time. Even with the cooling factor β = 0.9999, the wrong solutions were not fully eliminated in the case of the F _H objective function. In the case of F _H , cooling factor β = 0.9997 appeared to be enough to eliminate wrong solutions.

The hybrid algorithm fulfilled the expectations, delivering the correct solutions. Due to the different distribution of the applied objective functions in the search space, different parameters of the SA component appeared to be necessary for the proper identification of the constant reference points. In the case of F _K , a slower cooling (β = 0.9995) was necessary. In the case of the F _H objective function, a faster cooling of SA (β = 0.999) could be applied to get the proper solution.

Although the computation cost is not a crucial element of the assessment of analysed algorithms, a comparison of the processing time for all algorithms and for both objective functions gives valuable information about the effectiveness of particular identification procedures. Figure 6 Shows the approximate processing times for 1000-time analysis for three algorithms.

Figure 6:

Comparison of processing time for 1000-time analysis.

Another issue that should be discussed in the context of the identification approach proposed in this study is the probability of detection of the particular groups of stable points. When appropriate parameters of objective functions are selected, particular groups of stable points are associated with the corresponding minima in the search space, as can be seen in Figure 6. The values of the objective function in these minima will depend on the number of points creating the group and the level of their internal congruency. The most desired result of the identification procedure is selecting all the groups with possibly balanced probability. Excluding the HJ algorithm, it can be easily seen that the SA algorithm significantly prefers the dominating solution (group 1, 2, 3, 4, 5) due to a lower objective function. It arises from the features of the SA algorithm, which is rather dedicated to finding the global minimum. In the case of a hybrid algorithm, both groups are detected with a more balanced probability. As mentioned before, such a result is more desirable due to the expert’s deciding role in the final decision.

4 Summary and conclusions

Proper identification of the stable reference points in a monitoring network is a complex problem. It is especially concerning in cases when the monitoring network is large and the identification result is equivocal. The presented research is a continuation of the research described in [23], and it utilizes the idea of the solution space search and application of optimization algorithms. The newly proposed objective function formulation using robust weight functions appeared to be an efficient way to detect the subgroups of congruent point groups. Two variants of the objective function based on Kadaj’s and Huber’s robust weight functions were applied. Both of them showed sufficient efficiency in detecting subgroups of congruent reference points. The F _K objective function, based on Kadaj’s weight function, although more numerically costly, in some cases delivers a cleaner solution. The definite conclusion needs, however, extended research on more numerous and diversified examples.

It also appeared that hybridization of the algorithm was the right solution. The hybrid algorithm, combining the HJ and SA procedures similarly to SA, selected only the proper groups of stable points, but it was significantly faster. Another advantage of the hybrid algorithm is higher sensitivity for minority subgroups, which should also be considered as a potential final reference system for deformations.

Although the stability identification result was limited only to the assumed groups of stable points, it can be expected that in real cases, the identification result will not always be clear and obvious, even with the application of the SA or Hybrid algorithm. For this reason, a filtering procedure must be applied for the processing of the raw result. The development of the optimal filtering algorithm needs further research.