The combined processing of geomagnetic intensity vector projections and absolute magnitude measurements

A method for the combined processing measurements of the projections and the absolute magnitudes of geomagnetic field intensity vectors, based on mathematical technology of local approximation models, is proposed. The approach realized in this paper, based on the proposed method, provides an increase in quality of the measurements of the projections of the geomagnetic intensity vector. An algorithm for the two-stage combined processing of the measurements of projections and absolute magnitudes of geomagnetic field intensity vectors is developed. The operation of the combined processing algorithm was tested on models and observatory measurements. The estimates of the combined processing algorithm errors were obtained using statistical modelling. The reduction of the root-mean-square error values was achieved for the estimates of the projections of geomagnetic field intensity vectors.

1 Introduction

In this article, a method and an algorithm for combined processing measurements of the components (projections and the absolute magnitudes) of geomagnetic field intensity vectors is proposed. The approach realized in the paper, based on the proposed method, provides an increase in quality of the measurements of the projections of the geomagnetic intensity vector. The considered measurements are carried out by INTERMAGNET observatories equipped with systems of vector and scalar magnetometers; the definitive type data are used, containing systematic errors equal to zero (Mandea and Korte, 2011; INTERMAGNET, 2018). The measurement errors of vector and scalar magnetometers are represented by random, normally distributed errors with zero expectation and a predetermined variance. As usual, the measurement errors of the projections of geomagnetic field intensity vectors are significantly larger than the ones of the absolute vector magnitudes performed by the mentioned measurement devices. The formulation of the problem of reducing the noise root-mean-square (RMS) error values in the geomagnetic field intensity projection measurements is due to combined processing of the values of all its components.

In the following research, the following steps are outlined:

1. A method for combined processing of the measurements of projections and absolute magnitudes of geomagnetic field intensity vectors is formulated, based on formation of the sequences of piecewise-constant local models followed by their weighted averaging;

2. A two-stage algorithm for combined processing of the measurements of projections and absolute magnitudes of geomagnetic field intensity vectors is developed;

3. Testing of the algorithm on model and observatory data is performed;

4. The estimates of the algorithm errors, calculated using statistical modelling, are presented; the reduction of the RMS noise errors for the estimates of the geomagnetic field vector projections is proved.

The material in this research paper is intended for specialists (magnetologists) engaged in digital processing of geomagnetic field measurements. The need to reduce noise errors in estimates of the projections of the geomagnetic field intensity vectors measured by vector magnetometers arises in a number of technical and scientific applications. For instance, technogenic disturbances can affect the geomagnetic observatory hardware, affecting vector and scalar magnetometers differently: as a rule, the noise errors from vector magnetometers occurring due to such interference are greater than the noise errors from scalar magnetometers. The decrease in noise errors from vector magnetometers is necessary, e.g. for calculation of the gradients of the projections of the geomagnetic field intensity vectors in the navigation problems.

Nowadays the reduction of errors for vector magnetometers (with certain assumptions) is achieved using optimization of the calibration from scalar magnetometers (Merrayo et al., 2000; Olsen et al., 2003) or refinement of calibration characteristics (Soborov et al., 2008) based on special mathematical processing. In the measurement systems considered in this paper, in fact, parallel measurements are performed; a possible algorithm providing the decrease in errors for such measurements can be formed based on Kalman filters (Shakhtarin, 2008). However, due to the peculiarities of this problem, the construction of the resulting non-linear filters is associated with certain problems due to the inaccuracies of linearization and the accepted hypothesis concerning the type of initial intensity vector function. In other research (Soloviev et al., 2018), joint processing of vector and scalar magnetometer measurements aimed at improving the calibration accuracy of the so-called baseline, which only indirectly provides the considered reduction in errors. The combined processing of measurements of projections and absolute magnitudes of geomagnetic field intensity vectors proposed in this paper is significantly free from the mentioned problems.

2 A method for combined processing of measurements of projections and absolute magnitudes of geomagnetic field intensity vectors

Let H₁(Ti), H₂(Ti), H₃(Ti), and $H_{\mathrm{0}}\left(Ti\right)=\sqrt{H_{\mathrm{1}}^{\mathrm{2}}\left(Ti\right)+H_{\mathrm{2}}^{\mathrm{2}}\left(Ti\right)+H_{\mathrm{3}}^{\mathrm{2}}\left(Ti\right)}$ be the initial functions for the projections and absolute magnitudes of the geomagnetic field intensity vectors; we assume that Y₁(Ti), Y₂(Ti), Y₃(Ti), and Y₀(Ti) are their values registered by vector and scalar magnetometers, $i=\mathrm{0},\mathrm{1},\mathrm{\dots },N_{\mathrm{f}}-\mathrm{1}$; the sampling interval T=1 s; 1 s measurements from INTERMAGNET observatories are analysed in this study. For $n=\mathrm{0},\mathrm{\dots },\mathrm{3}$, we represent the noise errors of the measurement values W_n(Ti) in the form of uncorrelated, normally distributed random values with zero mathematical expectations and some variances. Such a representation of errors is, to a large extent, valid for cases of large technogenic noises that can occur when geomagnetic measurements are carried out. We consider, with some assumptions, that the spectrum of random components for the functions of the geomagnetic field is concentrated almost entirely in the low-frequency domain, and the spectrum of random measurement errors is concentrated in the high-frequency domain.

We assume that the measurements, the initial functions, and the errors are related by linear additive dependences:

$$Y_n\left(Ti\right)=H_n\left(Ti\right)+W_n\left(Ti\right),n=\mathrm{0},\mathrm{\dots },\mathrm{3}.$$

Using the specified observation values Y_n(Ti), we demand the determination of the estimates $Y_{\mathrm{1}}^{\circ }\left(Ti\right)$, $Y_{\mathrm{2}}^{\circ }\left(Ti\right)$, and $Y_{\mathrm{3}}^{\circ }\left(Ti\right)$, where $i=\mathrm{0},\mathrm{1},\mathrm{\dots },N_{\mathrm{f}}-\mathrm{1}$, that would be close to the initial functions of the intensity vector projections. We perform the combined processing for the projections and magnitudes of the geomagnetic field intensity vectors in two stages.

At the first stage, on the main interval with the points $i=\mathrm{0},\mathrm{1},\mathrm{\dots },N_{\mathrm{f}}-\mathrm{1}$, we introduce the N-point sliding local intervals with limiting points $N_{{\mathrm{1}}_j}$, $N_{{\mathrm{2}}_j}$, and the sliding step N_d as well as the quantity of sliding intervals m₀.

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& N_{{\mathrm{1}}_j}=N_{\mathrm{d}}\left(j-\mathrm{1}\right),N_{{\mathrm{2}}_j}=N_{{\mathrm{1}}_j}+N-\mathrm{1},\\ & j=\mathrm{1},\mathrm{\dots },m_{\mathrm{0}}.\end{array}\end{array}$$

To simplify the considerations, we require the relations of multiplicity mN=N_f and N_dm_d=N; here, if m and m_d are integers, then $m_{\mathrm{0}}=m_{\mathrm{d}}(m-\mathrm{1})+\mathrm{1}$.

For a sliding interval with a number j, we define the model functions of a form $Y_{M_{\mathrm{1}}}\left(c_{{\mathrm{1}}_j},Ti\right)$, $Y_{M_{\mathrm{2}}}\left(c_{{\mathrm{2}}_j},Ti\right)$, and $Y_{M_{\mathrm{3}}}\left(c_{{\mathrm{3}}_j},Ti\right)$; here, $c_{n_j}$ is the vectors of model parameters, n=1, 2, 3. These model functions can be, in particular, polynomial, piecewise constant, piecewise linear, etc. The size of local intervals determines the approximation errors. At small N there will be large fluctuation errors, and at large N there will be large systematic approximation errors.

Based on the above-defined measured values, models, and the maximum likelihood method (Kramer, 1975), we define the local functional S(c_jY_j) which determines the measure of closeness for local measurements and models, similar to Getmanov (2013) as the sum of the four functionals:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& S_{\mathrm{0}}\left(c_{{\mathrm{1}}_j},c_{{\mathrm{2}}_j},c_{{\mathrm{3}}_j},Y_{{\mathrm{0}}_j}\right)=\\ & \sum _{i=N_{{\mathrm{1}}_j}}^{N_{{\mathrm{2}}_j}}{\left(Y_{{\mathrm{0}}_j}\left(Ti\right)-\sqrt{\sum _{r=\mathrm{1}}^{\mathrm{3}}{Y}_{M_r}^{\mathrm{2}}\left(c_{r_j},Ti\right)}\right)}^{\mathrm{2}}\\ & S_n\left(c_{n_j},Y_{n_j}\right)=\sum _{i=N_{{\mathrm{1}}_j}}^{N_{{\mathrm{2}}_j}}{\left(Y_{n_j}\left(Ti\right)-Y_{M_n}\left(c_{n_j},Ti\right)\right)}^{\mathrm{2}},\\ & n=\mathrm{1},\mathrm{2},\mathrm{3},\\ & S\left(c_j,Y_j\right)=S_{\mathrm{0}}\left(c_{{\mathrm{1}}_j},c_{{\mathrm{2}}_j},c_{{\mathrm{3}}_j},Y_{{\mathrm{0}}_j}\right)+\sum _{n=\mathrm{1}}^{\mathrm{3}}{S}_n\left(c_{n_j},Y_{n_j}\right).\end{array}\end{array}$$

Here, $c_j^T=(c_{{\mathrm{1}}_j}^T{c}_{{\mathrm{2}}_j}^T,c_{{\mathrm{3}}_j}^T)$ and $Y_j^T=(Y_{{\mathrm{1}}_j}^T{Y}_{{\mathrm{2}}_j}^T,Y_{{\mathrm{3}}_j}^T)$ are the parameter and value vectors related to the jth local interval. Taking into account the assumption of errors in measurements, we carry out the identification of the optimal estimates of the model parameters $c_j^{\circ }$ using the solutions of the sequence of optimization problems for local functionals:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}& c_j^{\circ }=\arg \left\{S(c_j,Y_j)\right\},j=\mathrm{1},\mathrm{\dots },m_{\mathrm{0}},\\ & c^{\circ T}=(c_{\mathrm{1}}^{\circ T},c_{\mathrm{2}}^{\circ T},c_{\mathrm{3}}^{\circ T}).\end{array}\end{array}$$

We perform the construction of sliding local models in an obvious way by assuming n=1, 2, 3, and $j=\mathrm{1},\mathrm{\dots },m_{\mathrm{0}}$.

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& Y_{n_j}^{\circ }\left(Ti\right)=Y_{M_n}\left(c_{n_j}^{\circ },Ti\right),N_{{\mathrm{1}}_j}\le i\le N_{{\mathrm{2}}_j},\\ & Y_{n_j}^{\circ }\left(Ti\right)=\mathrm{0},\mathrm{0}\le i<N_{{\mathrm{1}}_j},N_{{\mathrm{2}}_j}<i\le N_{\mathrm{f}}-\mathrm{1}.\end{array}\end{array}$$

At the second stage we introduce the unit functions E_j(Ti), equal to zero outside of local sliding intervals, add them together to get E(Ti), and calculate the sequence of weighting coefficients R(Ti):

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& E_j\left(Ti\right)=\mathrm{1},N_{{\mathrm{1}}_j}\le i\le N_{{\mathrm{2}}_j},E_j\left(Ti\right)=\mathrm{0},\\ & \mathrm{0}\le i<N_{{\mathrm{1}}_j},N_{{\mathrm{2}}_j}<i\le N_{\mathrm{f}}-\mathrm{1};\\ & E\left(Ti\right)=\sum _{j=\mathrm{1}}^{m_{\mathrm{0}}}{E}_j\left(Ti\right),R\left(Ti\right)=\mathrm{1}/E\left(Ti\right),\\ & i=\mathrm{0},\mathrm{1},\mathrm{\dots },N_{\mathrm{f}}-\mathrm{1}.\end{array}\end{array}$$

Let us perform the weighting averaging using Eq. (5) for the sum of the sliding local model sequence (Getmanov et al., 2015).

$$\begin{array}{}\text{(6)}& Y_n^{\circ }\left(Ti\right)=R\left(Ti\right)\sum _{j=\mathrm{1}}^{m_{\mathrm{0}}}{Y}_{n_j}^{\circ }\left(Ti\right)\end{array}$$

We calculate the estimates $Y_n^{\circ }\left(Ti\right)$ and $i=\mathrm{0},\mathrm{1},\mathrm{\dots },N_{\mathrm{f}}-\mathrm{1}$, where n=1, 2, 3, based on linear operations in Eq. (6) corresponding to the second stage.

The method of combined two-step processing of the values of the projections and the absolute magnitudes of geomagnetic field intensity vectors consists of the sequential execution of the first and second stages in accordance with Eqs. (2)–(6).

3 An algorithm for two-stage processing of measurements of projections and absolute magnitudes of geomagnetic field intensity vectors for piecewise-constant models

Let us build the local intervals using Eq. (1) and define the local models on them as piecewise-constant functions $Y_{M_n}\left(c_{n_j},Ti\right)=c_{n_j}$, $j=\mathrm{1},\mathrm{\dots },m_{\mathrm{0}}$, and $N_{{\mathrm{1}}_j}\le i\le N_{{\mathrm{2}}_j}$. In this case it is obvious that the initial functions for geomagnetic field vector projections must be approximately constant at local intervals with duration NT. The local interval can be expanded if we treat piecewise-constant functions as local models. Let us formulate the equation for local functionals:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}S\left(c_j,Y_j\right)=& \sum _{n=\mathrm{1}}^{\mathrm{3}}\left(\sum _{i=N_{{\mathrm{1}}_j}}^{N_{{\mathrm{2}}_j}}{\left(Y_{n_j}\left(Ti\right)-c_{n_j}\right)}^{\mathrm{2}}\right.\\ & \left.+\sum _{i=N_{{\mathrm{1}}_j}}^{N_{{\mathrm{2}}_j}}{\left(Y_{{\mathrm{0}}_j}\left(Ti\right)-\sqrt{c_{{\mathrm{1}}_j}^{\mathrm{2}}+c_{{\mathrm{2}}_j}^{\mathrm{2}}+c_{{\mathrm{3}}_j}^{\mathrm{2}}}\right)}^{\mathrm{2}}\right).\end{array}\end{array}$$

We differentiate Eq. (7) with respect to $c_{n_j}$, equate the derivatives to zero, and get the necessary conditions for an extremum in the form of a system of three non-linear algebraic equations as a result:

$$\begin{array}{rl}{\displaystyle \frac{\partial S\left(c_j,Y_j\right)}{\partial c_{n_j}}}=& {\sum }_{i=N_{{\mathrm{1}}_j}}^{N_{{\mathrm{2}}_j}}{Y}_{n_j}\left(Ti\right)-\mathrm{2}N{c}_{n_j}\\ & +{\displaystyle \frac{c_{n_j}{\sum }_{i=N_{{\mathrm{1}}_j}}^{N_{{\mathrm{2}}_j}}{Y}_{{\mathrm{0}}_j}\left(Ti\right)}{\sqrt{c_{{\mathrm{1}}_j}^{\mathrm{2}}+c_{{\mathrm{2}}_j}^{\mathrm{2}}+c_{{\mathrm{3}}_j}^{\mathrm{2}}}}}=\mathrm{0},\\ & n=\mathrm{1},\mathrm{2},\mathrm{3}.\end{array}$$

In this case, an exact analytical solution to this system is possible. Omitting the calculations, we obtain expressions for local estimates, $j=\mathrm{1},\mathrm{\dots },m_{\mathrm{0}}$:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& c_{n_j}^{\circ }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}N}}\\ & \left(\mathrm{1}+{\displaystyle \frac{{\sum }_{r=N_{{\mathrm{1}}_j}}^{N_{{\mathrm{2}}_j}}{Y}_{{\mathrm{0}}_j}\left(Tr\right)}{\sqrt{{\left({\sum }_{r=N_{{\mathrm{1}}_j}}^{N_{{\mathrm{2}}_j}}{Y}_{{\mathrm{1}}_j}\left(Tr\right)\right)}^{\mathrm{2}}+{\left({\sum }_{r=N_{{\mathrm{1}}_j}}^{N_{{\mathrm{2}}_j}}{Y}_{{\mathrm{2}}_j}\left(Tr\right)\right)}^{\mathrm{2}}+{\left({\sum }_{r=N_{{\mathrm{1}}_j}}^{N_{{\mathrm{2}}_j}}{Y}_{{\mathrm{3}}_j}\left(Tr\right)\right)}^{\mathrm{2}}}}}\right)\\ & {\sum }_{i=N_{{\mathrm{1}}_j}}^{N_{{\mathrm{2}}_j}}{Y}_{n_j}\left(Ti\right),\end{array}\end{array}$$

where n=1, 2, 3, and get the piecewise-constant functions for local estimates $Y_n^{\circ }\left(Ti\right)=c_{n_j}^{\circ }$, $N_{{\mathrm{1}}_j}\le i\le N_{{\mathrm{2}}_j}$, $Y_n^{\circ }\left(Ti\right)=\mathrm{0}$, $\mathrm{0}\le i<N_{{\mathrm{1}}_j}$, $N_{{\mathrm{2}}_j}<i\le N_{\mathrm{f}}-\mathrm{1}$, and n=1, 2, 3, using Eq. (8) according to Eq. (4); let us present them as the realization of the first stage.

Weighted averaging of sequences of piecewise-constant local estimates and the calculations of the estimate functions $Y_n^{\circ }\left(Ti\right)$ for the second stage are performed using Eqs. (5) and (6).

4 Testing of the algorithm for combined processing on model and observatory measurements

4.1 Testing on model data

The developed algorithm for combined processing was tested on model measurements. Initial model functions for the projections of the geomagnetic intensity vector $H_{M{G}_n}\left(Ti\right)$ were presented as quadratic functions; the model function for the absolute vector magnitude $H_{M{G}_{\mathrm{0}}}\left(Ti\right)$ was calculated based on $H_{M{G}_n}\left(Ti\right)$, n=1, 2, 3.

$$\begin{array}{rl}& H_{M{G}_n}\left(Ti\right)=a_{\mathrm{0},n}+a_{\mathrm{1},n}Ti+a_{\mathrm{2},n}{\left(Ti\right)}^{\mathrm{2}},\\ & H_{M{G}_{\mathrm{0}}}\left(Ti\right)=\sqrt{H_{M{G}_{\mathrm{1}}}^{\mathrm{2}}\left(Ti\right)+H_{M{G}_{\mathrm{2}}}^{\mathrm{2}}\left(Ti\right)+H_{M{G}_{\mathrm{3}}}^{\mathrm{2}}\left(Ti\right)},\end{array}$$

and $i=\mathrm{0},\mathrm{1},\mathrm{\dots },N_{\mathrm{f}}-\mathrm{1}$. The noise error functions W_n(Ti) and n=0, 1, 2, 3 were modelled using normally distributed values with zero mathematical expectation and the variances ${\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}={\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}={\mathit{\sigma }}_{\mathrm{3}}^{\mathrm{2}}={\mathit{\sigma }}^{\mathrm{2}}$ and ${\mathit{\sigma }}_{\mathrm{0}}^{\mathrm{2}}$. The model values of the projections the absolute magnitude of the intensity vector were represented as sums:

$$\begin{array}{}\text{(9)}& \begin{array}{rl}& H_{M_n}\left(Ti\right)=H_{M{G}_n}\left(Ti\right)+W_n\left(Ti\right),\\ & i=\mathrm{0},\mathrm{1},\mathrm{\dots },N_{\mathrm{f}}-\mathrm{1},\end{array}\end{array}$$

where n=0, 1, 2, 3. The following parameter values for the model functions were assigned: $a_{\mathrm{0},\mathrm{1}}=\mathrm{30}\mathrm{000}$, $a_{\mathrm{1},\mathrm{1}}=\mathrm{0.4206}$, $a_{\mathrm{2},\mathrm{1}}=-\mathrm{0.0040}$, $a_{\mathrm{0},\mathrm{2}}=\mathrm{4300.4}$, $a_{\mathrm{1},\mathrm{2}}=a_{\mathrm{1},\mathrm{2}}=-\mathrm{0.4252}$, ${\mathit{\alpha }}_{\mathrm{2},\mathrm{2}}=\mathrm{0.0043}$, $a_{\mathrm{0},\mathrm{3}}=\mathrm{4089.5}$, $a_{\mathrm{1},\mathrm{3}}=\mathrm{0.5013}$, ${\mathit{\alpha }}_{\mathrm{2},\mathrm{3}}=-\mathrm{0.0056}$, and T=1 s. Model measurement values $H_{M_n}\left(Ti\right)$ were the input of Eq. (8) of the algorithm for combined processing based on piecewise-constant models.

The case of local intervals without sliding was considered, with the number of points N, $N_{\mathrm{1}j}=N(j-\mathrm{1})$ and $N_{\mathrm{2}j}=N_{\mathrm{1}j}+N-\mathrm{1}$, where $j=\mathrm{1},\mathrm{\dots },m$ and mN=N_f. The local estimates $H_{n_j}^{\circ }=c_j^{\circ }$, where $j=\mathrm{1},\mathrm{\dots },m$, were calculated, and the sequences of piecewise constant estimates $H_{L_n}^{\circ }\left(Ti\right)$, n=1, 2, 3, corresponding to the first processing stage, were formed.

For a sliding local intervals case with the number of points N, the sliding step N_d was selected, as well as the number of sliding local intervals m₀. Local estimates $H_{n_j}^{\circ }=c_j^{\circ }$, $j=\mathrm{1},\mathrm{\dots },m_{\mathrm{0}}$, and the estimates of functions $H_{n_j}^{\circ }\left(Ti\right)$, where n=1, 2, 3 and $i=\mathrm{0},\mathrm{1},\mathrm{\dots },N_{\mathrm{f}}-\mathrm{1}$, corresponding to the first processing stage, were calculated. Next, the second stage was performed where the estimates $H_{S_n}^{\circ }\left(Ti\right)$ were found.

Figure 1 Results of testing the processing algorithm on model measurement data: 1 – initial model functions, 2 – noised model functions, 3 – first-stage result (piecewise-constant model estimates without sliding), and 4 – second-stage result (estimates with weighted averaging).

For testing, the following values were selected: N_f=96, N=12, m=8, N_d=1, m₀=84, σ=1.0 nT, and σ₀=0.5 nT. In Fig. 1, the calculation results are displayed for $H_{M{G}_{\mathrm{1}}}$ (Fig. 2a) and $H_{M{G}_{\mathrm{2}}}$ (Fig. 2b), $H_{M{G}_{\mathrm{3}}}$ (Fig. 2c); dashed lines with index 1 depict the initial functions for the projections of the geomagnetic field intensity vector $H_{M{G}_n}\left(Ti\right)$; lines with index 2 represent the noised measurements of the projections of the geomagnetic field intensity $H_{M_n}\left(Ti\right)$; piecewise-constant lines with index 3 represent the results of the first stage $H_{L_n}^{\circ }\left(Ti\right)$; solid lines with index 4 show the estimates $H_{S_n}^{\circ }\left(Ti\right)$ – the second-stage results with weighted averaging.

Figure 2 Results of testing the processing algorithm on observatory measurement data: 1 – initial data, 2 – first-stage result ( piecewise constant model estimates without sliding), and 3 – second-stage result (estimates with weighted averaging with sliding).

The performed testing of the processing algorithm on model data for a number of parameters led to the conclusion that the second stage of processing reduces the RMS of the errors of the first stage by 60 %–80 % on average.

4.2 Testing on real observatory geomagnetic data

The developed algorithm was tested using combined processing of 1 s sampled geomagnetic measurements from the INTERMAGNET observatory MBO (Mbour, Senegal). The measurements were recorded on 2 January 2014, they began at 01:16:37 UT, and the length of a test fragment was 96 s (N_f=96). For the processing algorithm, N=12 and the sliding step N_d=1 were assigned.

In Fig. 2 the test results are shown for H₁ (Fig. 2a), H₂ (Fig. 2b), and H₃ (Fig. 2c). Dashed lines with index 1 depict the observatory measurements of the geomagnetic vector projections H_n(Ti); piecewise-constant lines with index 2 are related to the first processing stage – the functions for piecewise constant estimates $H_{L_n}^{\circ }\left(Ti\right)$ without sliding are displayed; index 3 stands for the line corresponding to the result of the second processing stage – the estimate with weighted averaging $H_{S_n}^{\circ }\left(Ti\right)$ with sliding.

The results of testing the algorithm for combined processing of measurements of projections and absolute magnitudes of geomagnetic field intensity, displayed on Figs. 1 and 2, proved their satisfactory performance.

5 Error estimation for the algorithm for combined processing of measurements

The estimates of errors of the proposed algorithm for combined processing were found using statistical modelling. The first stage of combined processing was analysed.

The initial functions for the intensity vector projections were assumed to be constant on a local interval. The values $H_n(k,m,H_{\mathrm{0}},Ti)=H_n(k,m,H_{\mathrm{0}})$, $i=\mathrm{0},\mathrm{1},\mathrm{\dots },N-\mathrm{1}$, n=1, 2, 3 were found using the following equations:

$$\begin{array}{rl}& H_{\mathrm{1}}(k,m,H_{\mathrm{0}})=H_{\mathrm{0}}\sin {\mathit{\vartheta }}_m,H_{\mathrm{2}}(k,m,H_{\mathrm{0}})=H_{\mathrm{0}}\cos {\mathit{\vartheta }}_m\sin {\mathit{\phi }}_k,\\ & H_{\mathrm{3}}(k,m,H_{\mathrm{0}})=H_{\mathrm{0}}\cos {\mathit{\vartheta }}_m\cos {\mathit{\phi }}_k.\end{array}$$

Here, H₀ is an assigned absolute magnitude value, φ_kϑ_m are the azimuthal and zenithal angles, θ_m=Δϑm, $\mathrm{\Delta }\mathit{\vartheta }=\mathrm{2}\mathit{\pi }/M_{\mathrm{0}}$, m and k are integer parameters, $m=\mathrm{0},\mathrm{1},\mathrm{\dots },M_{\mathrm{0}}-\mathrm{1}$, φ_k=Δφk, $\mathrm{\Delta }\mathit{\phi }=\mathrm{2}\mathit{\pi }/K_{\mathrm{0}}$, and $k=\mathrm{0},\mathrm{1},\mathrm{\dots },K_{\mathrm{0}}-\mathrm{1}$.

For all possible values of indices k and m, the realizations of sequences of model normally distributed random numbers with zero mathematical expectation $W_{n,s}(k,m,Ti)$ were formed, where $i=\mathrm{0},\mathrm{1},\mathrm{\dots },N-\mathrm{1}$, and n=0, 1, 2, 3, and $s=\mathrm{1},\mathrm{\dots },S_{\mathrm{0}}$ is the number of realization for statistical modelling. For n=0, the variance determining the noise error level for a scalar magnetometer assumed the value ${\mathit{\sigma }}_{\mathrm{0}}^{\mathrm{2}}$; for n=1, 2, 3, the noise error variances for a vector magnetometer assumed the value σ². For $H_n(k,m,H_{\mathrm{0}})$, $W_{n,s}(k,m,Ti)$, and S₀, random realizations were constructed:

$$\begin{array}{rl}& H_{n,s}(k,m,H_{\mathrm{0}},i)=H_n(k,m,H_{\mathrm{0}})+W_{n,s}(k,m,Ti),\\ & i=\mathrm{0},\mathrm{1},\mathrm{\dots },N-\mathrm{1},n=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{3},s=\mathrm{1},\mathrm{\dots },S_{\mathrm{0}}.\end{array}$$

The results of the algorithm implementation for the first stage – the estimates $H_{n,s}^{\circ }(k,m,H_{\mathrm{0}},N,\mathit{\sigma },{\mathit{\sigma }}_{\mathrm{0}})$, n=1, 2, 3, $k=\mathrm{0},\mathrm{1},\mathrm{\dots },K_{\mathrm{0}}-\mathrm{1}$, $m=\mathrm{0},\mathrm{1},\mathrm{\dots },M_{\mathrm{0}}-\mathrm{1}$, $s=\mathrm{1},\mathrm{\dots },and{S}_{\mathrm{0}}$ – were calculated depending on the parameters $H_{\mathrm{0}},N,\mathit{\sigma },{\mathit{\sigma }}_{\mathrm{0}}$. The error of the processing algorithm ${\mathit{\epsilon }}_n^{\mathrm{2}}(k,m,H_{\mathrm{0}},N,\mathit{\sigma },{\mathit{\sigma }}_{\mathrm{0}})$ was found using averaging over the number of realizations for fixed $n,k,$, and m:

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& {\mathit{\epsilon }}_n^{\mathrm{2}}(k,m,H_{\mathrm{0}},N,\mathit{\sigma },{\mathit{\sigma }}_{\mathrm{0}})={\displaystyle \frac{\mathrm{1}}{S_{\mathrm{0}}}}\sum _{s=\mathrm{1}}^{S_{\mathrm{0}}}\left(H_n(k,m,H_{\mathrm{0}})\right.\\ & {\left.-H_{n,s}^{\circ }(k,m,H_{\mathrm{0}},N,\mathit{\sigma },{\mathit{\sigma }}_{\mathrm{0}})\right)}^{\mathrm{2}}.\end{array}\end{array}$$

The error described by Eq. (10) was averaged over the number of the absolute geomagnetic vector projections n=1, 2, 3 and then over different k and m. The final formula for estimating the error was the following:

$$\begin{array}{}\text{(11)}& \begin{array}{rl}& {\mathit{\epsilon }}_{\mathrm{f}}^{\mathrm{2}}(H_{\mathrm{0}},N,\mathit{\sigma },{\mathit{\sigma }}_{\mathrm{0}})={\displaystyle \frac{\mathrm{1}}{K_{\mathrm{0}}{M}_{\mathrm{0}}}}\sum _{k=\mathrm{1}}^{K_{\mathrm{0}}}\sum _{m=\mathrm{1}}^{M_{\mathrm{0}}}{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\\ & \sum _{n=\mathrm{1}}^{\mathrm{3}}{\mathit{\epsilon }}_n^{\mathrm{2}}(k,m,H_{\mathrm{0}},N,\mathit{\sigma },{\mathit{\sigma }}_{\mathrm{0}}).\end{array}\end{array}$$

The results of the combined processing algorithm for the first stage were compared with the results of the operation of a possible linear filtering algorithm that was separately applied to the recordings of vector magnetometer channels. The linear filtering algorithm in this case was represented by standard equations:

$$\begin{array}{}\text{(12)}& H_{\mathrm{1}n,s}^{\circ }(k,m,H_{\mathrm{0}},N,\mathit{\sigma })={\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{0}}^{N-\mathrm{1}}{H}_{n,s}(k,m,H_{\mathrm{0}},i).\end{array}$$

The error estimate ${\mathit{\epsilon }}_{\mathrm{1}\mathrm{f}}^{\mathrm{2}}\left(H_{\mathrm{0}}N\mathit{\sigma }\right)$ for the linear filtering algorithm was calculated similar to Eqs. (10) and (11).

$$\begin{array}{}\text{(13)}& \begin{array}{rl}& {\mathit{\epsilon }}_{\mathrm{1}n}^{\mathrm{2}}(k,m,H_{\mathrm{0}},N,\mathit{\sigma })={\displaystyle \frac{\mathrm{1}}{S_{\mathrm{0}}}}\sum _{s=\mathrm{1}}^{S_{\mathrm{0}}}\left(H_n(k,m,H_{\mathrm{0}})\right.\\ & {\left.-H_{\mathrm{1}n,s}^{\circ }(k,m,H_{\mathrm{0}},N,\mathit{\sigma })\right)}^{\mathrm{2}}\\ & {\mathit{\epsilon }}_{\mathrm{1}\mathrm{f}}^{\mathrm{2}}\left(H_{\mathrm{0}}N\mathit{\sigma }\right)={\displaystyle \frac{\mathrm{1}}{K_{\mathrm{0}}{M}_{\mathrm{0}}}}\sum _{k=\mathrm{1}}^{K_{\mathrm{0}}}\sum _{m=\mathrm{1}}^{M_{\mathrm{0}}}{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\\ & \sum _{n=\mathrm{1}}^{\mathrm{3}}{\mathit{\epsilon }}_{\mathrm{1}n}^{\mathrm{2}}(k,m,H_{\mathrm{0}},N,\mathit{\sigma }).\end{array}\end{array}$$

The efficiency of the proposed algorithm for combined processing of measurements was estimated using the introduction of a relative decrease factor for the RMS error values $\mathit{\rho }(H_{\mathrm{0}},N,\mathit{\sigma },{\mathit{\sigma }}_{\mathrm{0}})$:

$$\begin{array}{}\text{(14)}& \mathit{\rho }(H_{\mathrm{0}},N\mathit{\sigma },{\mathit{\sigma }}_{\mathrm{0}})={\displaystyle \frac{{\mathit{\epsilon }}_{\mathrm{f}}^{\mathrm{2}}(H_{\mathrm{0}},N,\mathit{\sigma },{\mathit{\sigma }}_{\mathrm{0}})}{{\mathit{\epsilon }}_{\mathrm{1}\mathrm{f}}^{\mathrm{2}}(H_{\mathrm{0}},N,\mathit{\sigma })}}.\end{array}$$

For statistical modelling, the following values have been assumed: H₀=10 000 nT, N=5, M₀=50, K₀=50, S₀=100, σ₀=0.5, 0.3, 0.1, 0.03, and σ= 0.156, 0.312, 0.625, 1.25, 2.50, 5.00, 10.0. Figure 3 displays the results of the $\mathit{\rho }(H_{\mathrm{0}},N,\mathit{\sigma },{\mathit{\sigma }}_{\mathrm{0}})$ calculation depending on ${\log }_{\mathrm{2}}(\mathit{\sigma }/\stackrel{\mathrm{‾}}{\mathit{\sigma }})$ and $\stackrel{\mathrm{‾}}{\mathit{\sigma }}=\mathrm{0.156}$. Numbers 1, 2, 3, and 4 indicate the estimates for σ₀=0.5, 0.3, 0.1, and 0.03, respectively. The estimates of the introduced factor made it possible to get an idea of the effectiveness of the proposed processing.

Figure 3 Results of calculation of the relative decrease factors for the errors: 1 – σ₀=0.5; 2 – σ₀=0.3; 3 – σ₀=0.1; 4 – σ₀=0.03.

Analysis of the graphs shows that, for a fixed value σ, the introduced factor ρ decreases when σ₀ decreases, which is physically understandable. It is also seen that this factor tends to limit values with increasing σ. For a fixed H₀ , an increase in N leads to a decrease in the factor ρ. The performed statistical modelling for a wide range of parameters shows that the estimates for ρ are about 0.15–0.3, which indicates the efficiency of the proposed combined processing.

6 Conclusions

The proposed method for combined processing of the measurements of projections and absolute magnitudes of geomagnetic field intensity vectors and the corresponding two-stage algorithm developed appear to be satisfactorily workable. The testing of the developed combined processing algorithm on model and observatory measurement data proved its efficiency.

The approach realized in this paper, based on the proposed method, provides an increase in quality of the measurements of the projections of the geomagnetic intensity vector, allowing us to eliminate possible unwanted disturbances of artificial (anthropogenic) origin.

Statistical modelling for the developed algorithm for combined processing of measurements shows that at the first stage the relative decrease factor of RMS errors can reach values of 0.15–0.3, and at the second stage the decrease in the first-stage RMS errors can reach approximately 60 %–80 %.

Further reduction of the RMS noise errors can be implemented based on combined processing using local piecewise linear models for the values of projections and absolute magnitudes of geomagnetic field intensity vectors. The proposed combined processing algorithm can be implemented for many practically important tasks, in particular, when optimizing the operation of three-component accelerometer systems, three-component angular velocity sensors, or other three-component data arranged in a similar way. The technique allows for processing of the data with sampling intervals smaller than 1 s.