## 2 Problem Formulation
For this purpose, the classical geodetic adjustment method using the least squares method (GLSA) or the robust adjustment method (GRA), based on the so-called equivalent weighting matrix of measurements with gross error, can be used. The functional and stochastic models and the optimisation criteria for the two indicated geodetic adjustment methods could be formulated as follows [25; 26; 27]:  
\[\left. \begin{array}{c}\mathbf{V}_{\text{x}}(k+1)\text{=}\mathbf{A}(k+1) \widehat{\mathbf{d}}_{\text{x}}(k+1)+\mathbf{L}(k+1)\text{ --functional model}\\ \mathbf{C}_{\text{x}}(k+1)\text{=}\delta_{0}^{2}(k+1)\mathbf{P}_{\text{x}}(k +1)^{-1}\text{ --stochastic model}\\ \mathbf{V}_{\text{x}}^{\text{T}}(k+1)\mathbf{P}_{\text{x}}(k+1)\mathbf{V}_{ \text{x}}(k+1)\text{=}min\text{ --classic adjustment crit.}\\ \mathbf{V}_{\text{x}}^{\text{T}}(k+1)\widehat{\mathbf{P}}_{\text{x}}(k+1) \mathbf{V}_{\text{x}}(k+1)\text{=}min\text{ --robust adjustment crit.}\end{array} \right\} \tag{3}\]  
Given the assumptions (Equation (3)), after appropriate mathematical transformations, Relations 4 and 5 can be created for the estimation of the correction vector \(\widehat{\mathbf{d}}_{\text{x}}(k+1)\) of the predicted coordinates of the SV position (state vector \(\widehat{\mathbf{x}}(k+1)^{-}\), established a priori):  
\[\left. \begin{array}{c}\widehat{\mathbf{d}}_{x}(k+1)\text{=}-\left(\mathbf{ A}^{\text{T}}(\mathbf{k}+1)\mathbf{P}_{\text{x}}(\mathbf{k}+1)\mathbf{A}(\mathbf{k}+1) \right)^{-1}\\ \mathbf{A}^{\text{T}}(\mathbf{k}+1)\mathbf{P}_{\text{x}}(\mathbf{k}+1)\mathbf{ L}(k+1)-\text{ for the GLSA method,}\\ \widehat{\mathbf{d}}_{x}(k+1)\text{=}-\left(\mathbf{A}^{\text{T}}(\mathbf{k}+1) \widehat{\mathbf{P}}_{\text{x}}(\mathbf{k}+1)\mathbf{A}(\mathbf{k}+1) \right)^{-1}\\ \mathbf{A}^{\text{T}}(\mathbf{k}+1)\widehat{\mathbf{P}}_{\text{x}}(\mathbf{k}+1 )\mathbf{L}(k+1)\text{ -- for the GRA method,}\end{array} \right. \tag{5}\]  
where:  
\[\mathbf{A}(k+1)\text{=}\left[ \begin{array}{cc}\frac{\vartheta f_{D_{1}}}{ \vartheta x}&\frac{\vartheta f_{D_{1}}}{\vartheta y}\\ \frac{\vartheta f_{R_{1}}}{\vartheta x}&\frac{\vartheta f_{R_{1}}}{\vartheta y}\\ \frac{\vartheta f_{D_{2}}}{\vartheta x}&\frac{\vartheta f_{D_{2}}}{\vartheta y}\\ \frac{\vartheta f_{R_{2}}}{\vartheta x}&\frac{\vartheta f_{R_{2}}}{\vartheta y}\\ \vdots\\ \frac{\vartheta f_{D_{n}}}{\vartheta x}&\frac{\vartheta f_{D_{n}}}{\vartheta y}\\ \frac{\vartheta f_{R_{n}}}{\vartheta x}&\frac{\vartheta f_{R_{n}}}{\vartheta y} \end{array} \right]\text{--design matrix,}\]  
Let us now assume that, based on the results of the measurements \(COG(k)\) and \(SOG(k)\) obtained using SV's onboard navigation equipment, e.g., the Inertial Navigation System (INS) from VectorNav 'NS-100' [28] and the Doppler Velocity Log (DVL) from LinkQuest 'NavQuest 600' [29] (ignoring the effects of current and sea drift), at any time moment \((k)\), during the time interval \(\tau\), the DR method 1 is used to determine the position coordinate vector \(\mathbf{x}(k+1)\). Knowing the covariance matrix \(\mathbf{P}(k)\) of the initial (previous) position coordinate vector \(\mathbf{x}(k)\) and the covariance matrix \(\mathbf{P}(k+1/1)\) of the coordinate vector increment \(\mathbf{x}(k+1)-\mathbf{x}(k)\), it is possible to calculate the covariance matrix \(\mathbf{P}(k+1)\) of the coordinate vector of the final (next) position \(\mathbf{x}(k+1)\) according to the relation:  
\[\mathbf{P}(k+1)=\mathbf{P}(k)+\mathbf{P}(k+1/1)=[ \begin{array}{cc}\sigma_{x }^{2}&\sigma_{xy}\\ \sigma_{yx}&\sigma_{y}^{2}\end{array} ]. \tag{6}\]  
In turn, the covariance matrix obtained in this way, \(\mathbf{P}(k+1)\), can be used indirectly to assess the accuracy of the position vector determination, \(\mathbf{x}(k+1)\), carried out with the following:  
* Mean error \[M_{xy}=\sqrt{\sigma_{x}{}^{2}+\sigma_{y}{}^{2}},\] (7)  
Figure 2: Algorithm for state vector estimation \(\tilde{\mathbf{x}}(k+1)\) using classical geodetic adjustment.