Application of Square Msplit Estimation in Determination of Vessel Position in Coastal Shipping

Abstract The main aim of this paper is to assess the possibility of using non-conventional geodetic estimation methods in maritime navigation. The research subject of this paper concerns robust determination of vessel’s position using a method of parameters estimation in the split functional model (Msplit estimation). The studies performed will help in finding out if and in which situations the application of Msplit estimation as the method for determining vessel’s position is beneficial from the perspective of navigation safety. The results obtained were compared with the results of traditional estimation methods, i.e. least squares method and robust M-estimation.


INTRODUCTION
The second half of the 20th century and the beginning of the 21st century, despite financial troubles observed worldwide, are a period of a very dynamic development of global economy, including growth of maritime economy. It is visible in the actual constant increase of vessel movements in maritime ports and thus in the increased traffic of vessels in port approaches. Such state of affairs gave rise to an urgent need to improve navigation safety, especially on port approach courses as well as in limited basins (both in terms of their area and depth). Since the beginning of the 80s, vessel traffic management systems (VTS) aiming to ensure high safety level of navigation in basins covered by them have played an important role. Similarly to air traffic control, they influence vessels' traffic flow. In order to fulfil their objectives, it is crucial for the VTS station controllers to have an up-to-date and highly precise position of the observed vessel. To this end, the information on the current position is transmitted from the vessel to shore via e.g. an AIS system or radio waves directly by a watch officer. However, a coastal station controller has no information about the quality of the positioning performed, for example using satellite systems (see e.g., Specht et al. 2015, Specht andRudnicki 2016), or the location where the positioning system antenna is set. Therefore, in order to precise the vessel obtained coordinates on the sea it is possible to make own observations using coastal radar stations and determine the position of the observed object. VTS systems are usually equipped with several coastal radar stations that transmit more observations than are needed. In order to use these redundant observations effectively, the authors suggest using new non-conventional methods of observations adjustment applied in modern geodesy.
The significance of the subject of this research consists of the fact that it deals with vessels' changing position; it is impossible to repeat the measurement at the same position. Therefore, the selection of an optimal observation processing method is very important. Empirical analyses were conducted under the assumption that the observation vector y consists of observations obtained from the coastal radar stations. Therefore, it may be assumed that navigation observations are a realization of a random variable Y~P x . The research is particularly focused on a case when one observation y i has a probability distribution with a distribution parameter that is different from the remaining observations or belongs to an acceptable family of probability distributions but has a significantly different standard deviation. The source literature uses a term of "outliers" for measurements of such characteristics. Outliers may be of different nature, e.g. misreading a measurement result, temporary changes in the measurement environment parameters or improper calibration of measurement equipment. Such outliers are usually interpreted as gross measurement errors. It is generally known that in such a case it is impossible to use the least squares method to process measurement results (Guo et al. 2010). This is due to high sensitivity of the least squares method to outliers. Therefore, the effect of outliers on the final results of estimated values needs to be eliminated or at least minimized. Currently, the methods relating to robust M-estimation (Yang et al. 2002, Guo et al. 2010, Ge et al. 2013) are predominantly used to this end. Then, it is possible to limit the influence of outliers on the final adjustment with the use of various weight functions, such as Huber's, Hampel's or Danish. This kind of approach to the processing of observations from coastal radar stations was proposed and described in detail in research papers (Czaplewski 2004, Świerczyński, Czaplewski 2013, 2015. However, the use of calculation strategy referring to the principles of robust M-estimation is often associated with the risk of incorrect specification of control parameters. Therefore, there is a risk of incorrect identification of a random measurement error as a gross error. The a priori nature of control parameters may constitute a serious limitation of applicability of the method described above. From the theoretical perspective, occurrence of an observation contaminated by a gross error in the y vector is related to the occurrence of a realization of a random variable with the probability distribution different than Y~P x in the observation set. Therefore, it may be concluded that the observation vector may be a mixture of realizations of several random variables, whereas only one of them is acceptable. Such an assumption constitutes the basis for writing down a probabilistic model of gross errors in the following manner (Huber 1981, Yang et al. 2002, Wiśniewski 2009): where 0;1 τ ∈< >. In the expression (1), (1) P X is the probability distribution of measurement error with acceptable values, belonging to the distribution family P={P X(1) :X (1) ∈Θ}, where Θ is the space of parameters. However, realizations of inacceptable random variables are treated as observations that have the distribution ( 2) P X and belong to the distribution family P={P X(2) :X (2) ∈Θ}. Thus, the influence of observations having the probability distribution significantly different than (1) P X is limited in the process of establishing robust M-estimates.
Therefore, we can note that the problem of gross errors occurring in the observations may be solved by using the method of parameter estimation in a split functional model. Although the method of M split estimation is a generalization of a wide class of M-estimation methods, it may be analyzed in a context similar for the classic, robust M-estimation (Wiśniewski 2009, Ge et al. 2013). The parameter estimation method in the split functional model presumes that each observation y 1 may be a realization of one out of two competitive random variables (1) (1)Ỹ P X and ( 2) (2)Ỹ P X (Wiśniewski 2009, Zienkiewicz, Baryła 2015. As part of an experiment, it was assumed that "good" observations have a probability distribution of (1) P X . That is why estimation of competitive versions of a probability distribution parameter X should result in values of the vector (1) X free from outliers (observations with gross errors will be "assigned" to the vector (2) X ). As a result, the vector (1) X will contain reliable coordinates of a vessel.
Empirical analyses were conducted on the basis of simulated observations. The determined M split estimates were compared both with the results of the robust M-estimation and the classic method of the least squares. Since the theoretical grounds of the robust M-estimation and classic method of the least squares are greatly described in the literature, these descriptions are omitted further in this paper and the main focus is placed on the results of the numerical test at the end of the paper. The research conducted by the authors aims at determining whether the use of the M split estimation may help a controller at the VTS station to obtain more precise position coordinates of the observed vessel, which should significantly enhance navigation safety.

THEORETICAL FOUNDATIONS OF M SPLIT ESTIMATION
The estimation theory usually adopts the following functional observation model: y = F(X) -v = F(X) + ε (2) where: y∈ℜ n is an observation vector, F(X) is a vector function specifying the measured value, X∈ℜ n is a vector of determined parameters, ε∈ℜ n is a random vector of measurement errors and v∈ℜ n is a vector of theoretic observation corrections. It is also assumed that the measurement results are mutually independent, i.e.
, : cov( , ) 0 i j i j y y ∀ = ; therefore, the measurement results covariance matrix is a diagonal matrix with the following form: The assumption that each y 1 observation may be a realization of one out of two competitive random variables (1) (1)Ỹ P X and ( 2) (2)Ỹ P X , belonging to the family of probability distribution P = {P X(1) , P X(2) : X (1) , X (2) ∈ Θ} constitutes theoretical grounds for the M split estimation method. As a consequence of such an assumption, the classic functional model of geodetic observations (2) is split into two competitive models (Wiśniewski 2009) (1) (1) where: X (1) and X (2) -competitive versions of a parameter X, v (1) and v (1) -two versions of theoretical observation corrections vector concerning the same observation vector y. At the beginning, for convenience, it is assumed that the split is a natural linear functional model = − v AX y done into two competing models (1) (1) n,r is a known matrix of coefficients. It must be noted that a split of the classic functional model may be done for any number of functional models. Then, it is a case of an M split(q) estimation, where q is the number of splits of the classic functional model (2). Such generalization of the parameter estimation method in the split functional model was proposed and described in detail in the paper (Wiśniewski 2010).
Assigning an observation to a "suitable" functional model depends on an elementary split potential K in an observation y 1 (Wiśniewski 2009, 2010, Wiśniewski and Zienkiewicz 2016. The split potential may be interpreted as a certain measure determining the odds of assigning the observation to any available, competitive probability distribution (1) ( , P P P ∈ X X . The papers (Wiśniewski 2009(Wiśniewski , 2010 prove that defining an optimization problem based on the elementary split potential leads to searching of the minimum of the following aim function: where: It must be noted that the expression (5) may be treated as a generalization of the classic M-estimation aim function (cf. Huber 1964, Hampel et al. 1986). Therefore, the M split estimation is a particular type of an M-estimation method development. Now, one should assume that general components of the aim function (5) adopt a characteristic form for the method of the square M split estimation (Wiśniewski 2009, Zienkiewicz 2014) This particular version of the M split estimation aim function reflects a case when the adopted probabilistic models belong to a family of normal distributions (Wiśniewski 2009(Wiśniewski , 2010. The aim function minimum (7) is generally searched using the Newton method (Teunissen 1990, Wiśniewski 2009, 2010, Zienkiewicz 2014, Zienkiewicz and Baryła 2015). Determination of M split estimates may also be done with the gradient "zeroing". It should be noted that the function gradient (7) has to be determined with relation to two arguments of the aim function, i.e. for arguments (1) X and (2) X respectively. Then, it can be stated that M split estimates (1) X and (2) X resolving the following equation system are searched (1) (1) Therefore, by determining particular gradient forms for the aim function (7) the iteration process of determining square the M split estimates for 1,..., j m = may be written down in the following form In the expression (11) as (12) and (13) were written cross-weight matrices, where ( ) . It should be noted that adopting the aim function components as square functions, i.e.   (1) i v , respectively. As a consequence, occurrence of high value of (1) i v correction in the observation contaminated by gross error will result in the cross-weight function (13) "assigning" such an observation with the competitive functional model. Graphic interpretation of cross-weight function in square M split estimation is presented in Fig. 1. The process of determining the square M split estimates is of iteration nature due to the cross-weighting. The start points in the iteration process may be estimates obtained with the least squares method. The iteration process ends with obtaining such values of the square M split estimates (1) X and (2) X , for which In these points, the condition of existence of the aim function minimum, i.e. (1) (1) (2)( ; , ) = g y X X 0and (2) (1) (2)( ; , ) = g y X X 0 has also to be fulfilled.

SYSTEM OF OBSERVATION EQUATIONS FOR THE ANALYZED REASERCH PROBLEM
If the observed position is determined on the basis of a bearing, it can be easily noted than the vector function ( ) F X is not linear in the expressions (2) and (4) ( ) arctan Usually, for the sake of convenience, non-linear functions are brought to a linear form by expanding these functions into a Taylor series limited to first expressions of the expansion. Then, the linearized function ( ) F X may be written down in the following form where: X° -approximate value of parameter vector X , dX -increment vector for approximate coordinates. Therefore, components of the vector ( ) F X  are approximate measured values that are a vector function of approximate parameter vector X°. A linear version of the split functional model (4) may be written on this basis (16) where: dX and (2) dX -two versions of competitive increment vectors such that (1) For the so-defined systems of observation equations, the competitive vectors (1) dX may be determined based on the aim function (7) with function arguments (1) dX and (2) dX in the M split estimation methods (1) On the basis of the determined increment estimates (1) dX and (2) dX , the adjusted values of their respective parameters, i.e.ˆd = Coordinates of a vessel constitute parameters for the navigation measurement structures, including the measurement network (created with the use of coastal radar stations and bearings as observations) studied in this paper. Thus, the vectors , determined in the algorithm (17) will contain increments to approximate two-dimension coordinates of a vessel. It should be noted that only vector (1) dX will contain reliable information about vessel's position. This is due to an earlier premise that the measurement gross errors are "assigned" to the functional model (2) (2) d = + v A X L . Therefore, the estimate (2) dX will be contaminated with outliers.
In view of the foregoing, the matrix of the known factors A and the intercept vector L for the vector function (15) have general forms as follows: (19) where are approximate bearings values calculated on the basis of approximate parameter values. The procedure for determining the factor matrix A, the weight matrix P and the intercept vector L when processing radar observations was described in detail in the papers Czaplewski 2013, 2015).

NUMERICAL EXAMPLE
The presented numerical test relates to a real case where navigation observations are obtained from coastal radar stations working in the "VTS Gdańsk Bay" system. In this system, it is possible to make simultaneous bearings to a vessel from several coastal radar stations. Therefore, a vast scope of the data obtained with radiolocation equipment makes it possible to determine a precise position of a vessel with the use of estimation methods, and in particular the robust estimation method. Fig. 2 shows an example of the measurement structure used.

Fig. 2. Measurement structure using coastal radar stations
In the empirical analysis, the vessel's bearing was measured simultaneously from five coastal radar stations. Coordinates of these stations in the PL-UTM system are presented in Table  1. When presenting the results, the authors of this paper opted for a certain simplification and resigned from inserting units of measurement next to the analyzed values. Therefore, it must be noted that all horizontal coordinates X and Y are indicated in meters [m] and bearings NR in degrees [°]. In this example, a simulation of vessel's position [ , ] i i on the Gdańsk Bay was performed three times for research purposes. The bearing was simulated assuming a theoretical position of a vessel in three measurement periods, i.e. , and . For the simulation of random measurement errors, a random number generator in the MATLAB software was used with the assumption that observation errors have normal distribution ~(0, ) N σ for 0.5 σ = . It was assumed that the bearing medium error values were equal to the adopted standard deviation 0.5 i y m = . On the basis of 0.5 i y m = , the specific form of weight matrix was determined that was used for determining the LS and M split estimates as well as it initiated the determination of a robust M-estimate. The values of simulated bearings together with a theoretical vessel position are presented in Table 2. It was assumed in the experiment that the bearings made from coastal station Gdynia_KP Harbour Master were contaminated with a gross error of 10 e =  in each of three measurement periods (  Table 3 focuses on the procedure for determining robust estimates on the basis of vessel's position Z 1 . The robust M-estimates were determined on the basis of the following iteration process (for 1,..., j c = ) where: ( ) w v -weight matrix.
In this paper, individual elements ( ) The values l and g are interpreted as control parameters, which have a significant influence on the length and efficiency of the iteration process. Whereas = v v -initiating the determination process of the Danish estimate -are analyzed in the starting step of the robust M-estimation with the "original" weight matrix. It should be noted that, at the starting step of the iteration process, the standardized observation corrections determined for the second and third bearing clearly show that the observations NR 2 and NR 3 need to be treated with caution when determining vessel coordinates. For the calculations, t = 2.5, g = 2 and l = 1/t were assumed. Such control parameter values are the most frequently assumed in practical use of the Danish weight function (Wiśniewski 2014). As a result, the first iteration of the process (20) using the weight function (21) determined new, equivalent weights of these observations, i.e. and . It can be easily noted that, except for a bearing contaminated with a gross error, a correct observation NR 3 was also suppressed". This is due to a significant value of the adopted parameter l. It is a case of the so-called "hard" suppression of observation, where weights of the observations which had standardized corrections significantly beyond the permissible range 2.5; 2.5 i v ∆ < − > were processed uncompromisingly. It should be noted that the iteration process was prematurely ended also for this reason. The first iteration of determining Danish M-estimate is de facto the last iteration step. For this reason, the control parameter values should be selected in experiments for the particular geometrical measurement structure (Nowel 2016). Nonetheless, in case of navigation networks, such analyses are impossible to perform due to a dynamic position change of the studied subject (vessel). The results of slightly milder suppression of outliers were presented in the paper (Świerczyński and Czaplewski 2015). However, it should be emphasized that, despite classifying the observation NR 3 as an outlier, the vessel position determined with the use of the Danish weight function has to be deemed correct.
As stated earlier, the least squares method estimates were also used as the starting elements for the process of determining M split estimates. Then 0 ( )= w v P . By analyzing  (2) (1) i W Danish dX (1) dX (2)  observations homogenous in terms of accuracy, it can be noted that the matrix P 2 does not influence the values of determined estimates (1) dX and (2) dX . Therefore, in the analyzed example, square of the observation weight matrix does not influence cross-weighting. Thus, Table 3 presents cross-weighting matrices as (2) W and (1) W . Since the values 0 (1) dX and 0 (1) v are already known, there is no obstacle for the values of the vectors (2) dX and (2) v to be calculated according to the rules of the M split estimation in the zero step of iteration (see formula (17)). After the first iteration step is performed, it can be easily seen that the vector 1 (1) dX becomes robust to the influence of gross error contaminating the bearing NR 2 . The contaminated observation is assigned to a competitive functional model (2) (2) means that the calculation algorithm, even after first step, marks the bearing NR 2 as "uncertain" and limits its influence onto the determined estimate (1) dX . The results of the last iteration confirm that the M split estimation method faultlessly identified the outlier. By analyzing the components of cross-weight matrices 2(2) 0.00 W = and 2(1) 105.91 W = it can be clearly stated that the contaminated bearing was substantially excluded from the process of determining the estimate (1) dX . The value of M split estimate free from gross error is approximate to the robust M-estimate determined with the use of the Danish weight function. Thus, positions of the vessel determined on the basis of these two estimates are comparable.
The estimates of vessel coordinates were determined according to the foregoing procedure for three variants of vessel position. For the purposes of comparison, the determination process of LS, Danish and M split estimates was conducted for a non-contaminated set y and considered a gross error in the bearing NR 1 . However, due to a wide range of performed empirical analyses, the study description is limited to presenting the final results of adjustment shown in Table 4. The determined LS estimates can be treated as a correctly determined vessel position for the observation sets not contaminated with outliers. Thus, the estimate determined with the least squares method may constitute a point of reference for other estimates. It can be easily noted that contamination of the bearing NR 2 with gross error leads to a situation when the determined LS estimate does not provide reliable information in neither of the three variants of vessel's position. However, the determined robust M-estimates and the M split estimates correctly find the gross error in the bearing measurement network and limit its influence on the determined position of the vessel. It should be emphasized that reliable information on vessel's location is provided also by the M split estimate if the vector does not contain outliers. Graphical interpretation of the obtained results is presented in Fig. 3.

Fig. 3. Position of a measurement unit determined on the basis of neutral and robust estimation methods
From the practical perspective, determination of influence of gross error e on the determined estimates of vector ˆL S dX is an interesting issue. The results of empirical analysis concerning this issue are presented in Fig. 4. The authors analyzed two variants of gross error occurrence: in bearing NR 2 (Fig. 4a) and NR 3 (Fig. 4b). The LS, M split and robust M-estimates for several variants of value 20, 20 e ∈ − were determined in the test. The obtained results clearly show that along with the increase of the gross error absolute value, contamination of the estimate in the least squares method increases significantly, thus contributing to incorrect assessment of vessel's position. The M split and the robust M-estimates efficiently limit the influence of an outlier on the determined position of a vessel. Moreover, it should be noted that the obtained M split estimates show higher repeatability than the determined M-estimates in the presented example. Nonetheless, this does not influence the positive evaluation of both methods in detecting bearings contaminated with gross error.

CONCLUSION
The source literature provides many arguments that creating bearing networks based on coastal radar stations supports accurate determination of a vessel position in areas of navigational risks or heavy traffic areas, thus increasing navigational safety (e.g. Czaplewski 2013, 2015). The bearing network related to the operation of the "VTS Gdańsk Bay" system and also analyzed in this paper is an example of such a network. In this positioning strategy, the precision of position determination depends highly on the adopted estimation method. The most natural estimation method for such measurement structures is the least squares method. The aforementioned papers prove that the robust M-estimation methods, used, for instance, in geodetic adjustment calculus, can be used for aligning of bearing networks. This results from the specific nature of the navigation observations, which can be contaminated with gross errors of different sources.
This paper proposes the use of parameter estimation method in the split functional model for determining a reliable vessel position.
On the basis of the weight function and influence function properties, the estimation methods can be classified as neutral, robust and weak (Kadaj 1988). It should be noted that, according to this classification, the least squares method is a method neutral for outliers, the Danish method belongs to the group of robust methods and the M split estimation belongs to the family of weak estimation methods. Generally, in the adjustment calculus theory, the weak estimation methods consist in determining an estimate favoring the outliers (Cellmer 2014). Nonetheless, the performed theoretical and empirical analyses concerning both methods show that the parameter estimation method in the split functional model may be analyzed in the same context as the robust M-estimation methods. Therefore, the obtained results confirm the initial hypothesis that the use of M split estimation method allows determination of a reliable position of a vessel, despite observations in the measurement structure that are contaminated with gross error. Similar properties were observed during empirical analysis for the robust M-estimation method, however, the M split estimates show higher repeatability of the determined vessel's coordinates.