Mathematical foundation of Capon's method for planetary magnetic field analysis

Minimum variance distortionless projection, the so-called Capon method, serves as a powerful and robust data analysis tool when working on various kinds of ill-posed inverse problems. The method has not only successfully been applied to multipoint wave and turbulence studies in the context of space plasma physics, but it is also currently being considered as a technique to perform the multipole expansion of planetary magnetic fields from a limited data set, such as Mercury's magnetic field analysis. The practical application and limits of the Capon method are discussed in a rigorous fashion by formulating its linear algebraic derivation in view of planetary magnetic field studies. Furthermore, the optimization of Capon's method by making use of diagonal loading is considered.

1 Introduction

Nonlinear and adaptive filter techniques have a wide range of applications in geophysical and space science studies to find the most likely parameter set describing the measurement data or to decompose the data into a set of signals and noise. Above all, the minimum variance distortionless projection introduced by Capon (1969) (hereafter, Capon's method) has successfully been applied to multipoint data analyses for waves, turbulence fields and current sheets (Motschmann et al., 1996; Glassmeier et al., 2001; Narita et al., 2003, 2013; Contantinescu et al., 2006; Plaschke et al., 2008). The strength of Capon's method lies in the fact that the method performs a robust data fitting even when the spatial sampling or data amount is limited in the measurement (e.g., successfully applied to four-spacecraft data, Motschmann et al., 1996). Capon's method is currently being considered for planetary magnetic field studies in which the data (i.e., magnetic field samples) are more limited (e.g., data sampled on single orbits) and has recently been applied to simulated Mercury magnetic field data in view of the BepiColombo mission (Toepfer et al., 2020).

From a theoretical point of view there are several origins for the derivation of the method. The first derivation of Capon's method (Capon, 1969), constructed for the analysis of seismic waves, is based on the estimation of frequency–wavenumber spectra. Later on, this approach was reformulated in terms of matrix algebra (Motschmann et al., 1996). In light of mathematical statistics, Capon's estimator can be regarded as a special case of the maximum likelihood estimator (Narita, 2019). In this work the linear algebraic formulation of the method (Motschmann et al., 1996) with specific attention to magnetic field analysis is extended and the application of diagonal loading is discussed to improve the quality of data analysis with more justified applications and limits.

2 Motivation of Capon's method

The analysis of planetary magnetic fields is of great interest and one of the main tasks in space science. Here we pay special attention to the analysis of Mercury's internal magnetic field, which is one of the primary goals of the BepiColombo mission (Benkhoff et al., 2010). The magnetometer onboard the Mercury Planetary Orbiter (MPO) (Glassmeier et al., 2010) measures the magnetic field vectors ${\mathit{b}}^i={\left(b_x^i,b_y^i,b_z^i\right)}^T\in {\mathbb{R}}^{\mathrm{3}}$ at N data points xⁱ, $i=\mathrm{1},\mathrm{\dots },N$ along the orbit in the vicinity of Mercury. The magnetic fields around Mercury are considered to be a composition or superposition of internal fields generated by the dynamo process as well as crustal and induced fields, which are mainly dominated by dipole and quadrupole fields and external fields generated by the currents flowing in the magnetosphere. For Mercury the external fields contribute a significant amount to the total magnetic field within the magnetosphere (Anderson et al., 2011), and therefore a robust method is required for separating the internal fields from the total measured field. Yet, each component has to be properly modeled and parameterized when decomposing the field.

For example, when only data in current-free regions are analyzed, the planetary magnetic field is non-rotational and can be parameterized via the Gauss representation (Gauß, 1839). Within these current-free regions, the internal magnetic field can be expressed as the gradient of a scalar potential Φ, which can be expanded into a set of basis functions. Because of

$$\begin{array}{}\text{(1)}& {\partial }_{\mathit{x}}\cdot \mathit{B}=\mathrm{0},\end{array}$$

the scalar potential Φ has to satisfy

$$\begin{array}{}\text{(2)}& {\partial }_{\mathit{x}}^{\mathrm{2}}\mathrm{\Phi }=\mathrm{0},\end{array}$$

where ${\partial }_{\mathit{x}}^{\mathrm{2}}$ is the Laplacian. Choosing planetary centered coordinates with radius $r^i\in \left[R_{\mathrm{M}},\mathrm{\infty }\right)$, azimuth angle ${\mathit{\lambda }}^i\in \left[\mathrm{0},\mathrm{2}\mathit{\pi }\right]$ and polar angle ${\mathit{\theta }}^i\in \left[\mathrm{0},\mathit{\pi }\right]$, the solution of Eq. (2) is given by spherical harmonics so that the potential results in

$$\begin{array}{}\text{(3)}& \begin{array}{rl}\mathrm{\Phi }\left({\mathit{x}}^i\right)=& R_{\mathrm{M}}\sum _{l=\mathrm{1}}^{\mathrm{2}}({\displaystyle \frac{R_{\mathrm{M}}}{r^i}}{)}^{l+\mathrm{1}}\sum _{m=\mathrm{0}}^l\\ & \left[g_l^m\cos \right(m{\mathit{\lambda }}^i)+h_l^m\sin (m{\mathit{\lambda }}^i\left)\right]P_l^m\left(\cos \right({\mathit{\theta }}^i\left)\right)\end{array}\end{array}$$

for the planetary dipole and quadrupole fields. Here, R_M indicates the radius of Mercury and $P_l^m$ represents the Schmidt-normalized associated Legendre polynomials of degree l and order m. Within this series expansion a set of expansion coefficients $g_l^m$ and $h_l^m$ occurs, called internal Gauss coefficients. By constructing the Gauss coefficients in a vectorial fashion and defining the true coefficient vector as $\mathit{g}={\left(g_{\mathrm{1}}^{\mathrm{0}},g_{\mathrm{1}}^{\mathrm{1}},h_{\mathrm{1}}^{\mathrm{1}},g_{\mathrm{2}}^{\mathrm{0}},g_{\mathrm{2}}^{\mathrm{1}},h_{\mathrm{2}}^{\mathrm{1}},g_{\mathrm{2}}^{\mathrm{2}},h_{\mathrm{2}}^{\mathrm{2}}\right)}^T$, the magnetic field vectors at every data point xⁱ and the expansion coefficients are related via

$$\begin{array}{}\text{(4)}& {\displaystyle }{\mathit{b}}^i& {\displaystyle }=-{\partial }_{{\mathit{x}}^i}\mathrm{\Phi }\left({\mathit{x}}^i\right)+{\mathit{v}}^i+{\mathit{n}}^i\text{(5)}& {\displaystyle }& {\displaystyle }={\mathbf{H}}^i\mathit{g}+{\mathit{v}}^i+{\mathit{n}}^i,\end{array}$$

where the terms of the series expansion are arranged in the matrix ${\mathbf{H}}^i\left(r^i,{\mathit{\theta }}^i,{\mathit{\lambda }}^i\right)$, called the shape matrix. The vector vⁱ describes the parts that are not parameterized by the underlying model, e.g., the external parts, covering the range between the parameterized field and the measurement noise of the sensors, which is symbolized by the vector nⁱ. The measurement noise is neither correlated with the parameterized part Hⁱg nor with the non-parameterized part vⁱ.

Summarizing the magnetic field measurements for all N data points into a vector $\mathit{B}={\left({\left({\mathit{b}}^{\mathrm{1}}\right)}^T\mathrm{\dots }{\left({\mathit{b}}^N\right)}^T\right)}^T\in {\mathbb{R}}^{\mathrm{3}N}$, the field can be written as

$$\begin{array}{}\text{(6)}& \mathit{B}=\mathbf{H}\mathit{g}+\mathit{v}+\mathit{n},\end{array}$$

where $\mathit{n}={\left({\mathit{n}}^{\mathrm{1}}\mathrm{\dots }{\mathit{n}}^N\right)}^T\in {\mathbb{R}}^{\mathrm{3}N}$, $\mathit{v}={\left({\mathit{v}}^{\mathrm{1}}\mathrm{\dots }{\mathit{v}}^N\right)}^T\in {\mathbb{R}}^{\mathrm{3}N}$ and $\mathbf{H}={\left[{\mathbf{H}}^{\mathrm{1}}\mathrm{\dots }{\mathbf{H}}^N\right]}^T\in {\mathbb{R}}^{\mathrm{3}N\times G}$. G indicates the number of expansion coefficients.

Within Eq. (6) the magnetic field vector B and the shape matrix H, given by the underlying model, are known. The coefficient vector g is to be determined by data fitting. Since in most applications the number of known magnetic data points is much larger than the number of wanted expansion coefficients (G≪3N), H is a rectangular matrix in general. Furthermore, the non-parameterized parts of the field and the noise are unknown. Therefore, the direct inversion of Eq. (6) is impossible and g has to be estimated. In this case, Capon's method establishes a robust and useful tool to find the estimated solution for the expansion coefficients in Eq. (6).

Since the method does not require the orthogonality of the basis functions, it has a wider range of applications when decomposing the measured data into a set of superposed signals, especially when the number of data points is limited. For example, when the magnetic field data are measured on a dense grid in the vicinity of the planet, the Gauss coefficients can be estimated via integration of the data. But in the case of a limited data set those integrals cannot be evaluated.

3 Derivation of Capon's method

The following derivation of Capon's method is based on the linear algebraic formulation (Motschmann et al., 1996), which was formerly applied to the analysis of plasma waves in the terrestrial magnetosphere. Now we are focusing on the analysis of planetary magnetic fields.

As illustrated in the previous section, the magnetic field bⁱ measured at data point xⁱ in the vicinity of Mercury and the wanted expansion coefficients g are related via

$$\begin{array}{}\text{(7)}& {\mathit{b}}^i={\mathbf{H}}^i\mathit{g}+{\mathit{v}}^i+{\mathit{n}}^i\end{array}$$

or in a compact form for all N data points

$$\begin{array}{}\text{(8)}& \mathit{B}=\mathbf{H}\mathit{g}+\mathit{v}+\mathit{n},\end{array}$$

where the shape matrix H describes the underlying model.

For every data point xⁱ the noise vector nⁱ is assumed to be Gaussian with variance σ_n and zero mean so that 〈nⁱ〉=0 and $\langle {\mathit{n}}^{\mathit{i}}\circ {\mathit{n}}^{\mathit{i}}\rangle ={\mathit{\sigma }}_n^{\mathrm{2}}\mathbf{I}$, where I is the identity matrix. The angular brackets indicate averaging over an ensemble, e.g., different samples, realizations or measurements. Therefore, $\langle {\mathit{b}}^i\rangle ={\mathbf{H}}^i\mathit{g}+\langle {\mathit{v}}^i\rangle$ holds and equivalently

$$\begin{array}{}\text{(9)}& \langle \mathit{B}\rangle =\mathbf{H}\mathit{g}+\langle \mathit{v}\rangle ,\end{array}$$

since the model H and the true coefficient vector g are not affected by the averaging.

Because H is not always a square matrix but is in general a rectangular matrix with different sizes between rows and columns, the direct inversion of Eq. (9) is not guaranteed. Let us ignore for the moment the nonexistence of H⁻¹ and write down the equation

$$\begin{array}{}\text{(10)}& \mathit{g}={\mathbf{H}}^{-\mathrm{1}}\left(\langle \mathit{B}\rangle -\langle \mathit{v}\rangle \right).\end{array}$$

Despite its simplicity it is obviously incorrect. As H⁻¹ does not exist, let us look for another matrix w called the filter matrix, which follows the structure of this equation and in principle fulfills the solution of Eq. (9) with respect to g. Furthermore, the non-parameterized parts 〈v〉 are unknown, and therefore it is desirable to truncate these parts by the filter matrix. Capon's method is just the procedure to construct the filter matrix w and to calculate or, more precisely, to estimate g. To do so, some helpful quantities are introduced. In order to distinguish between the true coefficient vector g and the estimated solution, in the following Capon's estimator will be symbolized by g_C.

For the inversion of Eq. (6) it is useful to rewrite the vectors B and g_C in terms of a matrix representation. Thus, the data covariance matrix M and the coefficient matrix P are introduced as follows:

$$\begin{array}{}\text{(11)}& {\displaystyle }\mathbf{M}=\langle \mathit{B}\circ \mathit{B}\rangle ={\displaystyle \frac{\mathrm{1}}{Q}}\sum _{\mathit{\alpha }=\mathrm{1}}^Q{\mathit{B}}^{\mathit{\alpha }}\circ {\mathit{B}}^{\mathit{\alpha }}\in {\mathbb{R}}^{\mathrm{3}N\times \mathrm{3}N},\text{(12)}& {\displaystyle }\mathbf{P}=\langle {\mathit{g}}_{\mathrm{C}}\circ {\mathit{g}}_{\mathrm{C}}\rangle \in {\mathbb{R}}^{G\times G}.\end{array}$$

Here, Q indicates the number of measurements at each data point xⁱ, for example the number of flybys at each data point, and the circle ∘ symbolizes the outer product, which is defined by

$$\begin{array}{}\text{(13)}& \mathit{x}\circ \mathit{y}=\mathbf{x}\cdot {\mathbf{y}}^{†}\in {\mathbb{R}}^{n\times m}\end{array}$$

for any pair of vectors x∈ℝⁿ, y∈ℝ^m. The dagger † indicates the Hermitian conjugate and the dot stands for the multiplication of the matrices $\mathbf{x}\in {\mathbb{R}}^{n\times \mathrm{1}}$ and ${\mathbf{y}}^{†}\in {\mathbb{R}}^{\mathrm{1}\times m}$. Therefore, the diagonal of the matrix P contains the quadratic averaged components of the wanted estimator. It is important to note here that the covariance matrix M must be statistically averaged over several measurements. Otherwise the matrix is singular with a vanishing determinant (see Appendix A) and the further analysis cannot be achieved.

If the non-invertibility of the matrix B^α∘B^α is neglected, for every measurement α=1, …, Q an estimator ${\mathit{g}}_{\mathrm{C}}^{\mathit{\alpha }}$ for the true coefficient vector g can be determined so that

$$\begin{array}{}\text{(14)}& {\mathit{B}}^{\mathit{\alpha }}=\mathbf{H}{\mathit{g}}_{\mathrm{C}}^{\mathit{\alpha }}+{\mathit{v}}^{\mathit{\alpha }}+{\mathit{n}}^{\mathit{\alpha }}\end{array}$$

is valid. Thereby, each estimator deviates from the true coefficient vector by an error vector ${\mathit{\epsilon }}^{\mathit{\alpha }}=\mathit{g}-{\mathit{g}}_{\mathrm{C}}^{\mathit{\alpha }}$. Note that because of the non-invertibility of the matrix B^α∘B^α the single estimator ${\mathit{g}}_{\mathrm{C}}^{\mathit{\alpha }}$ cannot be calculated. Since the invertibility is solely given by the averaging over Q measurements, only the averaged estimator

$$\begin{array}{}\text{(15)}& {\mathit{g}}_{\mathrm{C}}={\displaystyle \frac{\mathrm{1}}{Q}}\sum _{\mathit{\alpha }=\mathrm{1}}^Q{\mathit{g}}_{\mathrm{C}}^{\mathit{\alpha }},\end{array}$$

with its related error

$$\begin{array}{}\text{(16)}& \langle \mathit{\epsilon }\rangle ={\displaystyle \frac{\mathrm{1}}{Q}}\sum _{\mathit{\alpha }=\mathrm{1}}^Q{\mathit{\epsilon }}^{\mathit{\alpha }},\end{array}$$

is available.

In contrast to the estimator, the true coefficient vector is a theoretical given vector that is not affected by the averaging (g≡〈g〉, Eq. 9). This property directly links the estimator to the true coefficient vector, which can be rewritten as

$$\begin{array}{}\text{(17)}& \mathit{g}={\mathit{g}}_{\mathrm{C}}^{\mathit{\alpha }}+{\mathit{\epsilon }}^{\mathit{\alpha }}.\end{array}$$

Averaging over Q measurements and using g≡〈g〉 results in

$$\begin{array}{}\text{(18)}& \mathit{g}={\mathit{g}}_{\mathrm{C}}+\langle \mathit{\epsilon }\rangle \end{array}$$

and analogously for the second-order moments

$$\begin{array}{}\text{(19)}& \mathit{g}\circ \mathit{g}=\langle {\mathit{g}}_{\mathrm{C}}\circ {\mathit{g}}_{\mathrm{C}}\rangle +\langle {\mathit{g}}_{\mathrm{C}}\circ \mathit{\epsilon }\rangle +\langle \mathit{\epsilon }\circ {\mathit{g}}_{\mathrm{C}}\rangle +\langle \mathit{\epsilon }\circ \mathit{\epsilon }\rangle .\end{array}$$

In the limit of vanishing errors 〈ε〉→0 and $\langle \mathit{\epsilon }\circ \mathit{\epsilon }\rangle \to \mathrm{0}$, Capon's estimator converges to the true coefficient vector:

$$\begin{array}{}\text{(20)}& {\mathit{g}}_{\mathrm{C}}\to \mathit{g},\end{array}$$

and therefore

$$\begin{array}{}\text{(21)}& \langle {\mathit{g}}_{\mathrm{C}}\circ {\mathit{g}}_{\mathrm{C}}\rangle \to \mathit{g}\circ \mathit{g}.\end{array}$$

For the further evaluation of Eq. (6) the definition of the outer product (Eq. 13) is utilized. Matrix multiplication of Eq. (6) with its Hermitian adjoint and averaging yields

$$\begin{array}{}\text{(22)}& \langle \mathit{B}\cdot {\mathit{B}}^{†}\rangle =\langle \left(\mathbf{H}\mathit{g}+\mathit{v}+\mathit{n}\right)\cdot {\left(\mathbf{H}\mathit{g}+\mathit{v}+\mathit{n}\right)}^{†}\rangle \end{array}$$

and therefore

$$\begin{array}{}\text{(23)}& \langle \mathit{B}\circ \mathit{B}\rangle =\mathbf{H}\cdot \mathit{g}\circ \mathit{g}\cdot {\mathbf{H}}^{†}+\mathrm{2}\left(\mathbf{H}\mathit{g}\right)\circ \langle \mathit{v}\rangle +\langle \mathit{v}\circ \mathit{v}\rangle +\langle \mathit{n}\circ \mathit{n}\rangle ,\end{array}$$

assuming that n is Gaussian with variance σ_n and zero mean (〈n〉=0). Because H g_C and 〈v〉 have the same dimension, the outer product commutates. By means of the limit 〈ε〉→0 in Eq. (21) the unknown matrix g∘g and the true coefficient vector g can be approximated by Capon's estimator, resulting in

$$\begin{array}{}\text{(24)}& \begin{array}{rl}\langle \mathit{B}\circ \mathit{B}\rangle & =\mathbf{H}\cdot \langle {\mathit{g}}_{\mathrm{C}}\circ {\mathit{g}}_{\mathrm{C}}\rangle \cdot {\mathbf{H}}^{†}+\mathrm{2}\left(\mathbf{H}{\mathit{g}}_{\mathrm{C}}\right)\circ \langle \mathit{v}\rangle \\ & +\langle \mathit{v}\circ \mathit{v}\rangle +\langle \mathit{n}\circ \mathit{n}\rangle .\end{array}\end{array}$$

Taking into account that $\langle \mathit{n}\circ \mathit{n}\rangle ={\mathit{\sigma }}_n^{\mathrm{2}}\mathbf{I}$ and using the above defined abbreviations, Eq. (24) can be rewritten as

$$\begin{array}{}\text{(25)}& \mathbf{M}=\mathbf{H}\mathbf{P}{\mathbf{H}}^{†}+\mathrm{2}\left(\mathbf{H}{\mathit{g}}_{\mathrm{C}}\right)\circ \langle \mathit{v}\rangle +\langle \mathit{v}\circ \mathit{v}\rangle +{\mathit{\sigma }}_n^{\mathrm{2}}\mathbf{I}\end{array}$$

Since Eq. (25) cannot be directly solved for P, the goal is to find the best estimator g_C for g so that $\mathbf{P}=\langle {\mathit{g}}_{\mathrm{C}}\circ {\mathit{g}}_{\mathrm{C}}\rangle$ is obtained as an approximate solution of Eq. (25). Therefore, a filter matrix w is constructed that separates the parameterized field from the noise by projecting the measured data onto the parameter space

$$\begin{array}{}\text{(26)}& {\mathbf{w}}^{†}\langle \mathit{B}\rangle ={\mathit{g}}_{\mathrm{C}}\end{array}$$

and simultaneously truncates the non-parameterized parts, i.e.,

$$\begin{array}{}\text{(27)}& {\mathbf{w}}^{†}\langle \mathit{v}\rangle =\mathrm{0}.\end{array}$$

Applying the filter matrix to the average of the non-parameterized parts of the field in Eq. (6),

$$\begin{array}{}\text{(28)}& \begin{array}{rl}\mathrm{0}& ={\mathbf{w}}^{†}\langle \mathit{v}\rangle ={\mathbf{w}}^{†}\left(\langle \mathit{B}\rangle -\mathbf{H}{\mathit{g}}_{\mathrm{C}}\right)={\mathbf{w}}^{†}\langle \mathit{B}\rangle -{\mathbf{w}}^{†}\mathbf{H}{\mathit{g}}_{\mathrm{C}}\\ & ={\mathit{g}}_{\mathrm{C}}-{\mathbf{w}}^{†}\mathbf{H}{\mathit{g}}_{\mathrm{C}},\end{array}\end{array}$$

where the true coefficient vector has been replaced by Capon's estimator (Eq. 18) and taking into account that 〈n〉=0 leads to the distortionless constraint

$$\begin{array}{}\text{(29)}& {\mathbf{w}}^{†}\mathbf{H}=\mathbf{I}.\end{array}$$

This equation is one of the important constraints for the construction of the wanted filter matrix w, but it is not enough. Let us look for another criterion. The basic idea is that in Eq. (6) there may be contributions 〈v〉 in the data B that are not caused by the internal magnetic field, and thus these contributions are not modeled by H g. Although the filter matrix w is already designed to truncate these parts, i.e., ${\mathbf{w}}^{†}\langle \mathit{v}\rangle =\mathrm{0}$, their contributions to the data are unknown, and therefore the parts that have to be truncated by w are unknown.

Referring to Eq. (26), the average output power, which is defined as the sum of the quadratic averaged components of the estimator, can be rewritten as (Pillai , 1989)

$$\begin{array}{}\text{(30)}& \mathrm{tr}\mathbf{P}=\mathrm{tr}\langle {\mathit{g}}_{\mathrm{C}}\circ {\mathit{g}}_{\mathrm{C}}\rangle =\mathrm{tr}\left({\mathbf{w}}^{†}\langle \mathit{B}\circ \mathit{B}\rangle \mathbf{w}\right),\end{array}$$

where tr 〈g_C∘g_C〉 is the trace of the matrix 〈g_C∘g_C〉.

Using Eq. (30), the coefficient matrix P can be expressed by the weight w and the data covariance M as

$$\begin{array}{}\text{(31)}& \mathbf{P}=\langle {\mathit{g}}_{\mathrm{C}}\circ {\mathit{g}}_{\mathrm{C}}\rangle ={\mathbf{w}}^{†}\langle \mathit{B}\circ \mathit{B}\rangle \mathbf{w}={\mathbf{w}}^{†}\mathbf{M}\mathbf{w}.\end{array}$$

The amount of noise and the amount of the non-parametrized part are unknown. Thus, it is conservative and safe to assume that a large part of the data is influenced by the noise and the non-parameterized parts that have to be truncated by the matrix w and the underlying model keeps the minimal contribution to the data. This minimal contribution has to be extracted.

Therefore, the output power P=tr P has to be minimized with respect to w^†, subject to the distortionless constraint ${\mathbf{w}}^{†}\mathbf{H}=\mathbf{I}$ or equivalently ${\mathbf{H}}^{†}\mathbf{w}=\mathbf{I}$. Using the Lagrange multiplier method this minimization problem can be formulated as

$$\begin{array}{}\text{(32)}& \mathrm{minimize}\mathrm{tr}\left[{\mathbf{w}}^{†}\mathbf{M}\mathbf{w}+\mathbf{\Lambda }\left(\mathbf{I}-{\mathbf{H}}^{†}\mathbf{w}\right)\right],\end{array}$$

where Λ represents the related Lagrange multipliers and the minimum is taken with respect to w. Since the components w_ij and $w_{ij}^{†}$ of the matrix w and w^†, respectively, are independent of each other, Eq. (32) can be expanded as

$$\begin{array}{}\text{(33)}& \mathrm{minimize}P=\mathrm{tr}\left[{\mathbf{w}}^{†}\mathbf{M}\mathbf{w}+\mathbf{\Lambda }\left(\mathbf{I}-{\mathbf{H}}^{†}\mathbf{w}\right)+\left(\mathbf{I}-{\mathbf{w}}^{†}\mathbf{H}\right)\mathbf{\Gamma }\right]\end{array}$$

or equivalently

$$\begin{array}{}\text{(34)}& \begin{array}{rl}\mathrm{minimize}P& =w_{ij}^{†}{M}_{jk}{w}_{ki}+{\mathrm{\Lambda }}_{ii}-{\mathrm{\Lambda }}_{ij}{H}_{jk}^{†}{w}_{ki}\\ & +{\mathrm{\Gamma }}_{ii}-w_{ij}^{†}{H}_{jk}{\mathrm{\Gamma }}_{ki}\end{array}\end{array}$$

with related additional Lagrange multipliers Γ. Taking the derivatives with respect to w_ki and $w_{ij}^{†}$ results in

$$\begin{array}{}\text{(35)}& \mathrm{0}={\partial }_{w_{ki}}P=w_{ij}^{†}{M}_{jk}-{\mathrm{\Lambda }}_{ij}{H}_{jk}^{†},\end{array}$$

yielding

$$\begin{array}{}\text{(36)}& {\mathbf{w}}^{†}\mathbf{M}=\mathbf{\Lambda }{\mathbf{H}}^{†}\end{array}$$

and

$$\begin{array}{}\text{(37)}& \mathrm{0}={\partial }_{w_{ij}^{†}}P=M_{jk}{w}_{ki}-H_{jk}{\mathrm{\Gamma }}_{ki},\end{array}$$

resulting in

$$\begin{array}{}\text{(38)}& \mathbf{M}\mathbf{w}=\mathbf{H}\mathbf{\Gamma }.\end{array}$$

Multiplication of Eq. (36) with w from the right and multiplication of Eq. (38) with w^† from the left considering the distortionless constraint delivers

$$\begin{array}{}\text{(39)}& \mathbf{P}={\mathbf{w}}^{†}\mathbf{M}\mathbf{w}=\mathbf{\Gamma }=\mathbf{\Lambda },\end{array}$$

and therefore

$$\begin{array}{}\text{(40)}& P=\mathrm{tr}\mathbf{P}=\mathrm{tr}\mathbf{\Gamma }=\mathrm{tr}\mathbf{\Lambda }.\end{array}$$

Because

$$\begin{array}{}\text{(41)}& P=\mathrm{tr}\mathbf{P}=|\langle {\mathit{g}}_{\mathrm{C}}\rangle {|}^{\mathrm{2}}\end{array}$$

is a convex function, Λ and Γ realize the minimal output power.

Due to the ensemble averaging the matrix M is invertible and M⁻¹ exists. Multiplying Eq. (38) with ${\mathbf{H}}^{†}{\mathbf{M}}^{-\mathrm{1}}$ and again considering the distortionless constraint yields

$$\begin{array}{}\text{(42)}& \mathbf{P}=\mathbf{\Lambda }=\mathbf{\Gamma }={\left[{\mathbf{H}}^{†}{\mathbf{M}}^{-\mathrm{1}}\mathbf{H}\right]}^{-\mathrm{1}}.\end{array}$$

By means of Eq. (36) the filter matrix results in

$$\begin{array}{}\text{(43)}& {\mathbf{w}}^{†}=\mathbf{P}{\mathbf{H}}^{†}{\mathbf{M}}^{-\mathrm{1}},\end{array}$$

and therefore Capon's estimator is given by

$$\begin{array}{}\text{(44)}& {\mathit{g}}_{\mathrm{C}}={\left[{\mathbf{H}}^{†}{\mathbf{M}}^{-\mathrm{1}}\mathbf{H}\right]}^{-\mathrm{1}}{\mathbf{H}}^{†}{\mathbf{M}}^{-\mathrm{1}}\langle \mathit{B}\rangle .\end{array}$$

Regarding the expensive derivation, this compact formula for Capon's estimator is surprising.

4 Diagonal loading

The filter matrix w is the key parameter to distinguish between the parameterized and the non-parameterized parts of the field. Conferring to Eq. (25) the ratio of these parts defines the input signal–noise ratio (given by the data) as

$$\begin{array}{}\text{(45)}& {\mathrm{SNR}}_{\mathrm{i}}={\displaystyle \frac{\mathrm{tr}\left(\mathbf{H}\mathbf{P}{\mathbf{H}}^{†}\right)}{\mathrm{tr}\left(\mathrm{2}\left(\mathbf{H}{\mathit{g}}_{\mathrm{C}}\right)\circ \langle \mathit{v}\rangle +\langle \mathit{v}\circ \mathit{v}\rangle +{\mathit{\sigma }}_n^{\mathrm{2}}\mathbf{I}\right)}}.\end{array}$$

The filter matrix is applied to the disturbed data for estimating the output power that is related to Capon's estimator. Thus, the output signal–noise ratio (resulting from the filtering) can be expressed as (Haykin, 2014; Van Trees, 2002)

$$\begin{array}{}\text{(46)}& \begin{array}{rl}{\mathrm{SNR}}_{\mathrm{o}}& ={\displaystyle \frac{\mathrm{tr}\left({\mathbf{w}}^{†}\mathbf{H}\mathbf{P}{\mathbf{H}}^{†}\mathbf{w}\right)}{\mathrm{tr}\left(\mathrm{2}{\mathbf{w}}^{†}\left(\mathbf{H}{\mathit{g}}_{\mathrm{C}}\right)\circ \langle \mathit{v}\rangle \mathbf{w}+{\mathbf{w}}^{†}\langle \mathit{v}\circ \mathit{v}\rangle \mathbf{w}\right)+{\mathit{\sigma }}_n^{\mathrm{2}}\mathrm{tr}\left({\mathbf{w}}^{†}\mathbf{w}\right)}}\\ & ={\displaystyle \frac{\mathrm{tr}\left(\mathbf{P}\right)}{{\mathit{\sigma }}_n^{\mathrm{2}}\mathrm{tr}\left({\mathbf{w}}^{†}\mathbf{w}\right)}},\end{array}\end{array}$$

since the filter fulfills the distortionless constraint ${\mathbf{w}}^{†}\mathbf{H}=\mathbf{I}$ and truncates the non-parameterized parts of the field, i.e., ${\mathbf{w}}^{†}\langle \mathit{v}\rangle =\mathrm{0}$. The ratio of the output and the input signal–noise ratio is the so-called array gain (Van Trees, 2002),

$$\begin{array}{}\text{(47)}& \begin{array}{rl}{\displaystyle \frac{{\mathrm{SNR}}_{\mathrm{o}}}{{\mathrm{SNR}}_{\mathrm{i}}}}=& {\displaystyle \frac{\mathrm{1}}{\mathrm{tr}\left({\mathbf{w}}^{†}\mathbf{w}\right)}}{\displaystyle \frac{\mathrm{tr}\left(\mathbf{P}\right)}{{\mathit{\sigma }}_n^{\mathrm{2}}\mathrm{tr}\left(\mathbf{H}\mathbf{P}{\mathbf{H}}^{†}\right)}}\\ & \mathrm{tr}\left(\mathrm{2}\left(\mathbf{H}{\mathit{g}}_{\mathrm{C}}\right)\circ \langle \mathit{v}\rangle +\langle \mathit{v}\circ \mathit{v}\rangle +{\mathit{\sigma }}_n^{\mathrm{2}}\mathbf{I}\right)\\ \sim & {\displaystyle \frac{\mathrm{1}}{\mathrm{tr}\left({\mathbf{w}}^{†}\mathbf{w}\right)}},\end{array}\end{array}$$

which is controlled by tr(w^†w), since H, P and v are given by the data and the underlying model.

The input signal–noise ratio is given by the data and the underlying model, and therefore ${\mathrm{SNR}}_{\mathrm{i}}=\mathrm{const}.$ in Eq. (47). In contrast to the input signal–noise ratio, the output signal–noise ratio depends on the filtering. When tr(w^†w) is large, the output signal–noise ratio can decrease (SNR_o→0), resulting in signal elimination, and thus the performance of Capon's estimator degrades. To prevent the signal elimination and to improve the robustness of Capon's estimator it is desirable to restrict tr(w^†w) with an upper boundary $T_{\mathrm{0}}=\mathrm{const}.$ (Van Trees, 2002), which can be expressed by the additional quadratic constraint

$$\begin{array}{}\text{(48)}& \mathrm{tr}\left({\mathbf{w}}^{†}\mathbf{w}\right)=T_{\mathrm{0}}.\end{array}$$

For reasons of mathematical aesthetics, the constant T₀ is expressed as the trace of a matrix T. For example, one can choose

$$\begin{array}{}\text{(49)}& \mathbf{T}={\displaystyle \frac{T_{\mathrm{0}}}{G}}\mathbf{I}\in {\mathbb{R}}^{G\times G},\end{array}$$

where G again indicates the number of wanted expansion coefficients and $\mathbf{I}\in {\mathbb{R}}^{G\times G}$ is the identity matrix so that

$$\begin{array}{}\text{(50)}& \mathrm{tr}\left(\mathbf{T}\right)={\displaystyle \frac{T_{\mathrm{0}}}{G}}\mathrm{tr}\left(\mathbf{I}\right)=T_{\mathrm{0}}\end{array}$$

holds. Thus, Eq. (48) can be rewritten as

$$\begin{array}{}\text{(51)}& \mathrm{tr}\left({\mathbf{w}}^{†}\mathbf{w}-\mathbf{T}\right)=\mathrm{0}.\end{array}$$

It should be noted that the matrix T can be chosen arbitrarily as long as it is independent of w and w^†.

Referring to the previous section and considering the additional quadratic constraint, the filter matrix can be calculated by solving

$$\begin{array}{}\text{(52)}& \begin{array}{rl}& \mathrm{minimize}\mathrm{tr}[{\mathbf{w}}^{†}\mathbf{M}\mathbf{w}+{\mathit{\sigma }}_{\mathrm{d}}^{\mathrm{2}}\left({\mathbf{w}}^{†}\mathbf{w}-\mathbf{T}\right)\\ & +\mathbf{\Lambda }\left(\mathbf{I}-{\mathbf{H}}^{†}\mathbf{w}\right)+\left(\mathbf{I}-{\mathbf{w}}^{†}\mathbf{H}\right)\mathbf{\Gamma }]\end{array}\end{array}$$

with respect to w, where ${\mathit{\sigma }}_{\mathrm{d}}^{\mathrm{2}}$ is the related additional Lagrange multiplier. Carrying out the same procedure as described in Sect. 3, the constrained minimizer is given by

$$\begin{array}{}\text{(53)}& \mathbf{w}={\left(\mathbf{M}+{\mathit{\sigma }}_{\mathrm{d}}^{\mathrm{2}}\mathbf{I}\right)}^{-\mathrm{1}}\mathbf{H}{\left[{\mathbf{H}}^{†}{\left(\mathbf{M}+{\mathit{\sigma }}_{\mathrm{d}}^{\mathrm{2}}\mathbf{I}\right)}^{-\mathrm{1}}\mathbf{H}\right]}^{-\mathrm{1}}.\end{array}$$

The comparison of Eq. (43) with Eq. (53) shows that the quadratic constraint results in the addition of the constant value ${\mathit{\sigma }}_{\mathrm{d}}^{\mathrm{2}}$ to the diagonal of the data covariance matrix M, which is known as diagonal loading (Van Trees, 2002). Consequently, the filter matrix is designed for a higher Gaussian background noise than is actually present (Van Trees, 2002).

In Fig. 1 tr(w^†w) is displayed with respect to $\mathit{\sigma }=\sqrt{{\mathit{\sigma }}_{\mathrm{d}}^{\mathrm{2}}+{\mathit{\sigma }}_n^{\mathrm{2}}}$. For σ_d→0, tr(w^†w) is large, resulting in a small output signal–noise ratio, which can cause signal elimination. For increasing values of σ_d, tr(w^†w) decreases and for σ_d→∞ Capon's filter converges to the least square fit filter

$$\begin{array}{}\text{(54)}& {\mathbf{w}}^{†}\stackrel{{\mathit{\sigma }}_{\mathrm{d}}\to \mathrm{\infty }}{⟶}{\left[{\mathbf{H}}^{†}\mathbf{H}\right]}^{-\mathrm{1}}{\mathbf{H}}^{†}\end{array}$$

or equivalently

$$\begin{array}{}\text{(55)}& {\mathbf{w}}^{†}\mathbf{w}\stackrel{{\mathit{\sigma }}_{\mathrm{d}}\to \mathrm{\infty }}{⟶}{\left[{\mathbf{H}}^{†}\mathbf{H}\right]}^{-\mathrm{1}},\end{array}$$

which treats all data equally.

Figure 1 Sketch of the trace of the filter matrix with respect to $\mathit{\sigma }=\sqrt{{\mathit{\sigma }}_{\mathrm{d}}^{\mathrm{2}}+{\mathit{\sigma }}_n^{\mathrm{2}}}$. For increasing values of the diagonal loading parameter σ_d the output signal–noise ratio increases and Capon's filter converges to the least square fit filter if the diagonal loading parameter is large.

Since tr(w^†w) is a monotonically decreasing function, σ_d→∞ might be the best choice for the diagonal loading parameter. To check this expectation we have to take a look at the output power P_d for the synthetic increased noise, which can be estimated by replacing $\mathbf{M}\to \mathbf{M}+{\mathit{\sigma }}_{\mathrm{d}}^{\mathrm{2}}\mathbf{I}$ in Eq. (42), resulting in (Pajovic, 2019; Richmond et al., 2005)

$$\begin{array}{}\text{(56)}& \begin{array}{rl}{P}_{\mathrm{d}}& =\mathrm{tr}\left[{\mathbf{w}}^{†}\left(\mathbf{M}+{\mathit{\sigma }}_{\mathrm{d}}^{\mathrm{2}}\mathbf{I}\right)\mathbf{w}\right]\\ & =\mathrm{tr}{\left[{\mathbf{H}}^{†}{\left(\mathbf{M}+{\mathit{\sigma }}_{\mathrm{d}}^{\mathrm{2}}\mathbf{I}\right)}^{-\mathrm{1}}\mathbf{H}\right]}^{-\mathrm{1}}\stackrel{{\mathit{\sigma }}_{\mathrm{d}}\to \mathrm{\infty }}{⟶}{\mathit{\sigma }}_{\mathrm{d}}^{\mathrm{2}}\mathrm{tr}{\left[{\mathbf{H}}^{†}\mathbf{H}\right]}^{-\mathrm{1}}.\end{array}\end{array}$$

Thus, P_d is an increasing function of σ_d. Since the output power has to be minimized one would expect that σ_d=0 is the best choice, which is in contradiction to the argumentation above. Therefore, the maximization of the array gain is not equivalent to the minimization of the output power (Van Trees, 2002). Since σ_d cannot be directly calculated within the minimization procedure (Eq. 52), it is not clear how to choose the optimal diagonal loading parameter σ_opt. that lies somewhere between those extrema.

In the literature there are several methods for determining the optimal diagonal loading parameter (Pajovic, 2019). In contrast to the analysis of waves, we favor measurements that are stationary up to the Gaussian noise for the analysis of planetary magnetic fields. Comparing the measurement times with planetary geology timescales, this assumption is surely valid for the internal magnetic field. For the external parts of the field this can be realized by choosing data sets of preferred situations, e.g., calm solar wind conditions. By virtue of the stationarity, the data covariance matrix $\mathbf{M}=\langle \mathit{B}\rangle \circ \langle \mathit{B}\rangle +{\mathit{\sigma }}_n^{\mathrm{2}}\mathbf{I}$ contains only one nontrivial eigenvalue ${\mathit{\lambda }}_{\mathrm{1}}=|\langle \mathit{B}\rangle {|}^{\mathrm{2}}+{\mathit{\sigma }}_n^{\mathrm{2}}$ and ${\mathit{\lambda }}_{\mathrm{i}}={\mathit{\sigma }}_n^{\mathrm{2}}$ for $i=\mathrm{2},\mathrm{\dots },\mathrm{3}\mathrm{N}$ elsewhere (see Appendix B). Therefore, estimators for the diagonal loading parameter that are related to eigenvalues corresponding to interference and noise (Carlson, 1988) or estimators taking into account the standard deviation of the diagonal elements of the data covariance matrix, i.e., ${\mathit{\sigma }}_{\mathrm{d}}^{\mathrm{2}}=\mathrm{std}\left\{{\mathit{\lambda }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\lambda }}_{\mathrm{3}\mathrm{N}}\right\}$ (Ma and Goh, 2003), cannot be applied. When simulated data are analyzed the deviation between Capon's estimator and the true coefficient vector, implemented in the simulation, is a useful metric for estimating the optimal diagonal loading parameter (Pajovic, 2019; Toepfer et al., 2020). But when the method is applied to real spacecraft data no true coefficient vector is available, and therefore another estimation method for the diagonal loading parameter that solely depends on the data and the underlying model is required.

The additional quadratic constraint (Eq. 48), resulting in diagonal loading, bounds the trace of the filter matrix that is the solution of the minimization procedure (Eq. 52). To prevent signal elimination, tr(w^†w) and P_d have to be minimized simultaneously. Since tr(w^†w) is decreasing for higher values of σ_d (see Fig. 1), P_d is an increasing function, and thus they act as competitors. At the optimal value σ_opt. the two competitors compromise, which can be understood as tr(w^†w) reaches its minimal value under the constraint that P_d is minimal. Therefore, the diagonal loading of the data covariance matrix is equivalent to the Tikhonov regularization for ill-posed problems. The Tikhonov regularization improves the robustness of the least square fit problem,

$$\begin{array}{}\text{(57)}& \mathrm{minimize}|\mathbf{H}\mathit{g}-\mathit{B}{|}^{\mathrm{2}},\end{array}$$

with respect to g under the constraint of solutions g with minimal norm, i.e.,

$$\begin{array}{}\text{(58)}& \mathrm{minimize}|\mathbf{H}\mathit{g}-\mathit{B}{|}^{\mathrm{2}}+\mathit{\alpha }|\mathit{g}{|}^{\mathrm{2}},\end{array}$$

where α is the corresponding Lagrange multiplier, which is also known as the regularization parameter (Tikhonov et al., 1995). In analogy to the Lagrange multiplier α, the optimal diagonal loading parameter can be estimated by the method of the L curve (Hiemstra et al., 2002). The L curve arises by plotting lg [tr(w^†w)] versus $\lg \left[\mathrm{tr}\left({\mathbf{w}}^{†}\left(\mathbf{M}+{\mathit{\sigma }}_{\mathrm{d}}^{\mathrm{2}}\mathbf{I}\right)\mathbf{w}\right)\right]$ for different values of σ_d and is displayed in Fig. 2. The optimal value σ_opt. is located in the vicinity of the L curve's knee, which is defined by the maximum curvature of the L curve (Hiemstra et al., 2002).

Figure 2 Sketch of the L curve for estimating the optimal diagonal loading parameter σ_d. The optimal value is located in the vicinity of the L curve's knee (dashed circle).

5 Practical application of Capon's method

When Capon's method is applied to the reconstruction of Mercury's internal magnetic field (Toepfer et al., 2020), several computationally burdensome matrix inversions are necessary to calculate the estimator (see Eq. 44)

$$\begin{array}{}\text{(59)}& {\mathit{g}}_{\mathrm{C}}={\left[{\mathbf{H}}^{†}{\mathbf{M}}^{-\mathrm{1}}\mathbf{H}\right]}^{-\mathrm{1}}{\mathbf{H}}^{†}{\mathbf{M}}^{-\mathrm{1}}\langle \mathit{B}\rangle .\end{array}$$

Thus, for the practical application of Capon's method it is useful to rewrite the method in terms of the least square fit method.

Capon's estimator inherits the matrix operation structure of the least square fit estimator, which is given by Haykin (2014):

$$\begin{array}{}\text{(60)}& {\mathit{g}}_L={\left[{\mathbf{H}}^{†}\mathbf{H}\right]}^{-\mathrm{1}}{\mathbf{H}}^{†}\langle \mathit{B}\rangle .\end{array}$$

Substituting $\mathbf{H}\to {\mathbf{M}}^{-\mathrm{1}/\mathrm{2}}\mathbf{H}$ and $\langle \mathit{B}\rangle \to {\mathbf{M}}^{-\mathrm{1}/\mathrm{2}}\langle \mathit{B}\rangle$ in Eq. (60), the least square fit estimator may be converted into Capon's estimator.

The least square fit method minimizes the deviation between the disturbed measurements B and the underlying model H g, measured in the Euclidian norm $\left|\right|.|{|}_{\mathrm{2}}$, so that

$$\begin{array}{}\text{(61)}& {\mathit{g}}_L=\arg \underset{\mathit{g}}{\min }\left|\right|\mathbf{H}\mathit{g}-\mathit{B}|{|}_{\mathrm{2}}^{\mathrm{2}}.\end{array}$$

Referring to the abovementioned substitutions, Capon's method can be interpreted as measuring the deviation in Eq. (61) in a different norm,

$$\begin{array}{}\text{(62)}& \left|\right|\mathbf{H}\mathit{g}-\mathit{B}|{|}_{{\mathbf{M}}^{-\mathrm{1}}}^{\mathrm{2}}=|\left|{\mathbf{M}}^{-\mathrm{1}/\mathrm{2}}\left(\mathbf{H}\mathit{g}-\mathit{B}\right)\right|{|}_{\mathrm{2}}^{\mathrm{2}},\end{array}$$

so that Capon's estimator is given by

$$\begin{array}{}\text{(63)}& {\mathit{g}}_{\mathrm{C}}=\arg \underset{\mathit{g}}{\min }\left|\right|\mathbf{H}\mathit{g}-\mathit{B}|{|}_{{\mathbf{M}}^{-\mathrm{1}}}^{\mathrm{2}}.\end{array}$$

Thus, Capon's method can be regarded as a special case of the least square fit method or more precisely of the weighted least square fit method, whereby the data and the model are measured and weighted with the inverse data covariance matrix. This property is useful for the practical application of Capon's method. In contrast to the computationally burdensome matrix inversions in Eq. (44), the usage of minimization algorithms, such as the gradient descent or conjugate gradient method, for solving Eq. (63) is more stable and computationally inexpensive.

To illustrate the mathematical foundations presented above, Capon's method is applied to simulated magnetic field data in order to reconstruct Mercury's internal magnetic field. As a proof of concept, the underlying model is restricted to the internal dipole and quadrupole contributions to the magnetic field as discussed in Sect. 2. Since simulated magnetic field data are analyzed, the special application to Mercury's magnetic field is of minor importance because the ideal solution is known from the simulation. If the method were tested against in situ measurement data, the application to the analysis of the Earth's magnetic field would be more suitable since the Earth's magnetic field is better known than Mercury's magnetic field.

For the reconstruction of Mercury's internal dipole and quadrupole field the internal Gauss coefficients $g_{\mathrm{1}}^{\mathrm{0}}=-\mathrm{190}\mathrm{nT}$ and $g_{\mathrm{2}}^{\mathrm{0}}=-\mathrm{78}\mathrm{nT}$ (Anderson et al., 2012; Wardinski et al., 2019), defining the nonvanishing components of the true coefficient vector g, are implemented in the hybrid code AIKEF (Müller et al., 2011), and the magnetic field data resulting from the plasma interaction of Mercury with the solar wind are simulated. The data are evaluated along selected parts of the prospective trajectories of the BepiColombo mission on the night side of Mercury within a distance of 0.2 R_M up to 0.4 R_M from Mercury's surface. Since simulated data are analyzed, the deviation between the true coefficient vector g and Capon's estimator g_C can be used as a metric to verify the estimation of the optimal diagonal loading parameter by making use of the L-curve technique. The optimal diagonal loading parameter results in ${\mathit{\sigma }}_{\mathrm{opt}.}\approx \mathrm{800}\mathrm{nT}$, which corresponds to the vicinity of the L curve's knee. The reconstructed Gauss coefficients for the internal dipole and quadrupole field are presented in Table 1.

Table 1 Implemented and reconstructed Gauss coefficients for the internal dipole and quadrupole field.

The deviation $|{\mathit{g}}_{\mathrm{C}}-\mathit{g}|$ between Capon's estimator and the true coefficient vector results in 30.2  nT or 14.7 %, respectively, for the optimal diagonal loading parameter. When the magnetic field data are evaluated at an ensemble of data points at a distance of 0.2 R_M from Mercury's surface this deviation is of the same order (Toepfer et al., 2020). Since the underlying model neglects the external parts and only the internal parts are considered, the coefficients $g_{\mathrm{2}}^{\mathrm{1}}$ and $h_{\mathrm{2}}^{\mathrm{1}}$ show large deviations from the implemented coefficients. Extending the underlying model by a parameterization of the external parts of the magnetic field improves Capon's estimator, especially when the magnetic field data are evaluated at some distance above Mercury's surface. In principle, one can also analyze the external field if some model is adopted.

6 Discussion of Capon's method

Capon's method has been previously applied to the analysis of waves. Thus, the existing derivations treat the non-parameterized parts of the field v as Gaussian noise so that v and the measurement noise n are of the same character (Motschmann et al., 1996). Considering the analysis of planetary magnetic fields this assumption is indefensible. For example, the external parts of the field are systematic noise and cannot be modeled by a Gaussian distribution. When the non-parameterized parts are Gaussian, i.e., 〈v〉=0, the truncation of these parts by the filter matrix ${\mathbf{w}}^{†}\langle \mathit{v}\rangle =\mathrm{0}$ (see Eq. 27) is fulfilled trivially, which reduces the terms within the derivation and the estimation of the diagonal loading parameter. Therefore, the mathematical foundations presented above generalize the previous derivations of Capon's method and transist into the derivation presented by Motschmann et al. (1996) for the special case of 〈v〉=0.

As already mentioned in the Introduction (Sect. 1), Capon's method can be regarded from several mathematical perspectives. Within the derivation presented above, the output power P, which is defined as the trace of the coefficient (covariance) matrix, is minimized with respect to the filter matrix w, subject to the distortionless constraint ${\mathbf{w}}^{†}\mathbf{H}=\mathbf{I}$. This procedure corresponds to the name minimum variance distortionless response estimator (MVDR), since P contains the variance of the model coefficients. Narita (2019) showed that Capon's estimator can also be derived by treating Capon's method as a special case of the maximum likelihood estimator by regarding the likelihood function as nearly Gaussian (particularly around the peak of the likelihood function). As discussed in Sect. 5, Capon's method can also be interpreted as a special case of the weighted least square fit. This illustrates that the several existing inversion methods for linear inverse problems are connected to each other and are not as different as they seem to be at first appearance.

7 Conclusions

Capon's method is a robust and useful tool for various kinds of ill-posed inverse problems, such as Mercury's planetary magnetic field analysis. The derivation of the method can be regarded from different mathematical perspectives. Here we revisited the linear algebraic matrix formulation of the method and extended the derivation for Mercury's magnetic field analysis. Capon's method becomes even more robust by incorporating the diagonal loading technique. Thereby, the construction of a filter matrix is vital to the derivation of Capon's estimator.

The trace of the filter matrix in particular determines the array gain, which is defined as the ratio of the output and the input signal–noise ratio. If the trace is large, the output signal–noise ratio can decrease, resulting in signal elimination, and thus the performance of Capon's estimator degrades. Bounding the trace of the filter matrix results in diagonal loading of the data covariance matrix, which improves the robustness of the method. The main problem of the diagonal loading technique is that in general it is not clear how to choose the optimal diagonal loading parameter. Since for the analysis of planetary magnetic fields we prefer measurements that are stationary up to the Gaussian noise, estimators for the diagonal loading parameter that are related to eigenvalues of the data covariance matrix corresponding to interference and noise cannot be applied. Making use of the L-curve technique enables a robust procedure for estimating the optimal diagonal loading parameter.

For the calculation of Capon's estimator several computationally burdensome matrix inversions are necessary. Interpretation of Capon's method as a special case of the least square fit method enables the usage of numerically more stable and fewer burden minimization algorithms, e.g., the gradient descent or conjugate gradient method, for calculating Capon's estimator.

It should be noted that the parameterization of Mercury's internal magnetic field via the Gauss representation, as mentioned in Sect. 2, is only one of several possibilities for modeling the magnetic field in the vicinity of Mercury. The underlying model can be expanded by other parameterizations, for example the Mie representation (toroidal–poloidal decomposition) or magnetospheric models, and Capon's method can be applied to estimate the related model coefficients. Besides the analysis of Mercury's internal magnetic field, the extension of the model also enables the reconstruction of current systems flowing in the magnetosphere. Concerning the BepiColombo mission this work establishes a mathematical basis for the application of Capon's method to analyze Mercury's internal magnetic field in a robust and manageable way.

Appendix A: Determinant of the outer product

The influence of the averaging on the determinant of the data covariance matrix M is exemplarily illustrated for the three-dimensional case. Thus, the magnetic field vector is given by

$$\begin{array}{}\text{(A1)}& \mathit{B}=\left(\begin{array}{c}{B}_x\\ B_y\\ B_z\end{array}\right),\end{array}$$

and therefore the outer product results in

$$\begin{array}{}\text{(A2)}& \mathit{B}\circ \mathit{B}=\left[\begin{array}{ccc}{B}_x^{\mathrm{2}}& B_x{B}_y& B_x{B}_z\\ B_x{B}_y& B_y^{\mathrm{2}}& B_y{B}_z\\ B_x{B}_z& B_y{B}_z& B_z^{\mathrm{2}}\end{array}\right]\end{array}$$

with a vanishing determinant

$$\begin{array}{}\text{(A3)}& \det \left(\mathit{B}\circ \mathit{B}\right)=\mathrm{3}{B}_x^{\mathrm{2}}{B}_y^{\mathrm{2}}{B}_z^{\mathrm{2}}-\mathrm{3}{B}_x^{\mathrm{2}}{B}_y^{\mathrm{2}}{B}_z^{\mathrm{2}}=\mathrm{0}.\end{array}$$

Throughout the averaging of the data, the data covariance matrix results in

$$\begin{array}{ll}{\displaystyle }\mathbf{M}& {\displaystyle }=\langle \mathit{B}\circ \mathit{B}\rangle =\langle \mathit{B}\rangle \circ \langle \mathit{B}\rangle +{\mathit{\sigma }}_n^{\mathrm{2}}\mathbf{I}\\ \text{(A4)}& {\displaystyle }& {\displaystyle }=\left[\begin{array}{ccc}\langle B_x^{\mathrm{2}}\rangle +{\mathit{\sigma }}_n^{\mathrm{2}}& \langle B_x\rangle \langle B_y\rangle & \langle B_x\rangle \langle B_z\rangle \\ \langle B_x\rangle \langle B_y\rangle & \langle B_y^{\mathrm{2}}\rangle +{\mathit{\sigma }}_n^{\mathrm{2}}& \langle B_y\rangle \langle B_z\rangle \\ \langle B_x\rangle \langle B_z\rangle & \langle B_y\rangle \langle B_z\rangle & \langle B_z^{\mathrm{2}}\rangle +{\mathit{\sigma }}_n^{\mathrm{2}}\end{array}\right]\end{array}$$

with a nonvanishing determinant

$$\begin{array}{}\text{(A5)}& \begin{array}{rl}\det \left(\mathbf{M}\right)& =(\langle B_x^{\mathrm{2}}\rangle +{\mathit{\sigma }}_n^{\mathrm{2}})(\langle B_y^{\mathrm{2}}\rangle +{\mathit{\sigma }}_n^{\mathrm{2}})(\langle B_z^{\mathrm{2}}\rangle +{\mathit{\sigma }}_n^{\mathrm{2}})\\ & -\langle B_x^{\mathrm{2}}\rangle \langle B_y^{\mathrm{2}}\rangle \langle B_z^{\mathrm{2}}\rangle -{\mathit{\sigma }}_n^{\mathrm{2}}\langle B_y^{\mathrm{2}}\rangle \langle B_z^{\mathrm{2}}\rangle \\ & -{\mathit{\sigma }}_n^{\mathrm{2}}\langle B_x^{\mathrm{2}}\rangle \langle B_z^{\mathrm{2}}\rangle -{\mathit{\sigma }}_n^{\mathrm{2}}\langle B_x^{\mathrm{2}}\rangle \langle B_y^{\mathrm{2}}\rangle \\ & ={\mathit{\sigma }}_n^{\mathrm{4}}(\langle B_x^{\mathrm{2}}\rangle +\langle B_y^{\mathrm{2}}\rangle +\langle B_z^{\mathrm{2}}\rangle )+{\mathit{\sigma }}_n^{\mathrm{6}}\\ & \ne \mathrm{0}.\end{array}\end{array}$$

Thus, the inverse of M exists, whereas the outer product
B∘B is singular.

Appendix B: Eigenvalues of the data covariance matrix

The data covariance matrix is defined as

$$\begin{array}{}\text{(B1)}& \mathbf{M}=\langle \mathit{B}\circ \mathit{B}\rangle =\langle \mathit{B}\rangle \circ \langle \mathit{B}\rangle +{\mathit{\sigma }}_n^{\mathrm{2}}\mathbf{I}.\end{array}$$

This matrix is quadratic and especially diagonalizable. Thus, there is a matrix D_M that is similar to the matrix M so that

$$\begin{array}{}\text{(B2)}& {\mathbf{D}}_{\mathbf{M}}={\mathbf{V}}^{†}\mathbf{M}\mathbf{V},\end{array}$$

where V is an orthogonal transformation, i.e., ${\mathbf{V}}^{†}\mathbf{V}=\mathbf{I}$, and D_M is a diagonal matrix, the diagonal elements of which are given by the eigenvalues of M. Inserting the definition of the matrix M delivers

$$\begin{array}{}\text{(B3)}& \begin{array}{rl}{\mathbf{D}}_{\mathbf{M}}& ={\mathbf{V}}^{†}\mathbf{M}\mathbf{V}={\mathbf{V}}^{†}\langle \mathit{B}\circ \mathit{B}\rangle \mathbf{V}\\ & ={\mathbf{V}}^{†}\left(\langle \mathit{B}\rangle \circ \langle \mathit{B}\rangle \right)\mathbf{V}+{\mathit{\sigma }}_n^{\mathrm{2}}\mathbf{I},\end{array}\end{array}$$

since ${\mathbf{V}}^{†}\mathbf{V}=\mathbf{I}$. To determine the diagonal form of the outer product, the two-dimensional case in which

$$\begin{array}{}\text{(B4)}& \langle \mathit{B}\rangle =\left(\begin{array}{c}\langle B_x\rangle \\ \langle B_y\rangle \end{array}\right)\end{array}$$

is considered exemplarily. Thus,

$$\begin{array}{}\text{(B5)}& \langle \mathit{B}\rangle \circ \langle \mathit{B}\rangle =\left[\begin{array}{cc}\langle B_x{\rangle }^{\mathrm{2}}& \langle B_x\rangle \langle B_y\rangle \\ \langle B_x\rangle \langle B_y\rangle & \langle B_y{\rangle }^{\mathrm{2}},\end{array}\right]\end{array}$$

and the characteristic polynomial results in

$$\begin{array}{}\text{(B6)}& \left(\langle B_x{\rangle }^{\mathrm{2}}-\mathit{\beta }\right)\left(\langle B_y{\rangle }^{\mathrm{2}}-\mathit{\beta }\right)-\langle B_x{\rangle }^{\mathrm{2}}\langle B_y{\rangle }^{\mathrm{2}}=\mathrm{0}\end{array}$$

or equivalently

$$\begin{array}{}\text{(B7)}& {\mathit{\beta }}^{\mathrm{2}}-\mathit{\beta }\left(\langle B_x{\rangle }^{\mathrm{2}}+\langle B_y{\rangle }^{\mathrm{2}}\right)={\mathit{\beta }}^{\mathrm{2}}-\mathit{\beta }|\langle \mathit{B}\rangle {|}^{\mathrm{2}}=\mathrm{0},\end{array}$$

where β denotes the eigenvalue given by β=0 and $\mathit{\beta }=|\langle \mathit{B}\rangle {|}^{\mathrm{2}}$. Therefore, in general the diagonal matrix of the outer product is given by

$$\begin{array}{}\text{(B8)}& {\displaystyle }{\mathbf{D}}_{\langle \mathit{B}\rangle \circ \langle \mathit{B}\rangle }{\displaystyle }=\left[\begin{array}{cccc}|\langle \mathit{B}\rangle {|}^{\mathrm{2}}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{⋮}& & \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{0}& \mathrm{\cdots }& \mathrm{\cdots }& \mathrm{0}\end{array}\right].\end{array}$$

The noise matrix ${\mathit{\sigma }}_n^{\mathrm{2}}\mathbf{I}$ is already a diagonal matrix so that the diagonal form of M is given by

$$\begin{array}{ll}{\displaystyle }{\mathbf{D}}_{\mathbf{M}}& {\displaystyle }=\left[\begin{array}{cccc}|\langle \mathit{B}\rangle {|}^{\mathrm{2}}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{⋮}& & \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{0}& \mathrm{\cdots }& \mathrm{\cdots }& \mathrm{0}\end{array}\right]+{\mathit{\sigma }}_n^{\mathrm{2}}\left[\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{⋮}& & \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{0}& \mathrm{\cdots }& \mathrm{\cdots }& \mathrm{1}\end{array}\right]\\ \text{(B9)}& {\displaystyle }& {\displaystyle }=\left[\begin{array}{cccc}|\langle \mathit{B}\rangle {|}^{\mathrm{2}}+{\mathit{\sigma }}_n^{\mathrm{2}}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{0}& {\mathit{\sigma }}_n^{\mathrm{2}}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{⋮}& & \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{0}& \mathrm{\cdots }& \mathrm{\cdots }& {\mathit{\sigma }}_n^{\mathrm{2}}\end{array}\right].\end{array}$$

Thus, the data covariance matrix contains only one nontrivial eigenvalue ${\mathit{\lambda }}_{\mathrm{1}}=\langle \mathit{B}\rangle {|}^{\mathrm{2}}+{\mathit{\sigma }}_n^{\mathrm{2}}$ and ${\mathit{\lambda }}_{\mathrm{i}}={\mathit{\sigma }}_n^{\mathrm{2}}$ for $i=\mathrm{2},\mathrm{\dots },\mathrm{3}N$.