Skip to main content

BRIEF RESEARCH REPORT article

Front. Phys., 09 June 2021
Sec. Space Physics

Error Propagation of Capon’s Minimum Variance Estimator

S. Toepfer
S. Toepfer1*Y. Narita,Y. Narita2,3D. HeynerD. Heyner3U. Motschmann,U. Motschmann1,4
  • 1Institut für Theoretische Physik, Technische Universität Braunschweig, Braunschweig, Germany
  • 2Space Research Institute, Austrian Academy of Sciences, Graz, Austria
  • 3Institut für Geophysik und Extraterrestrische Physik, Technische Universität Braunschweig, Braunschweig, Germany
  • 4DLR Institute of Planetary Research, Berlin, Germany

The error propagation of Capon’s minimum variance estimator resulting from measurement errors and position errors is derived within a linear approximation. It turns out, that Capon’s estimator provides the same error propagation as the conventionally used least square fit method. The shape matrix which describes the location depence of the measurement positions is the key parameter for the error propagation, since the condition number of the shape matrix determines how the errors are amplified. Furthermore, the error resulting from a finite number of data samples is derived by regarding Capon’s estimator as a special case of the maximum likelihood estimator.

1 Introduction

The reconstruction of model parameters from a given set of measurements is one of the most important tasks in geophysical and space science studies. The measurements are always affected by measurement errors which result in estimation errors for the wanted model parameters. The Capon method [13], also known as minimum variance distortionless response estimator (MVDR), is currently being considered as a robust inversion method for the analysis of planetary magnetic fields. In the past, the method has successfully been applied to the analysis of waves [1, 4] and therefore, specific attention has been paid to errors of the spectrum resulting from random perturbations in the amplitude and phase of sensor arrays [5]. Since Capon’s method is based on the evaluation of statistically averaged data, the spectrum is also affected by errors resulting from a finite number of samples [6]. Within the estimation of the frequency-wavenumber spectrum, the difference between static structures and waves or their combination can be discerned through dispersion relation analysis [712]. Concerning the application of the Capon method for the analysis of planetary magnetic fields, the error propagation of Capon’s estimator itself is of major importance for assessing the quality of the reconstructed model parameters. In general, two essential types of errors are expectable. On the one hand, the measured magnetic field data are affected by e.g., offsets, gains resulting from thermal variations and spacecraft magnetic disturbances [13]. On the other hand, the determination of the spacecraft’s positions can be defective (measurement position errors), which results in a defective shape matrix. As a follow-up of the generalized derivation of Capon’s method [3] and the error estimation of the power spectrum [6], in this work the effects of measurement errors, measurement position errors as well as finite sample sizes on Capon’s estimator are considered.

2 Capon’s Method

Before deducing the error propagation of Capon’s method, the main ideas of the method are shortly revisited [2, 3]. Due to the complexity of several physical problems, the entire parametrization of experimental data is unrealizable. Thus, it is useful to decompose the measurements B¯ into parametrized parts H¯¯g¯, where g¯ contains the corresponding wanted model coefficients and the shape matrix H¯¯ describes the distribution of the measurement positions with respect to the underlying model, non-parametrized parts v¯ as well as measurement noise n¯, so that

B¯=H¯¯g¯+v¯+n¯(1)

is valid. The measurement noise is assumed to be Gaussian with variance σn and zero mean (n¯=0). Since the shape matrix H¯¯ is not invertible and the non-parametrized parts are unknown, the exact solution for the wanted model coefficients g¯ is not available in general. Capon’s method delivers an estimator g¯C for the ideal solution g¯. The method is based on the construction of a filter matrix w¯¯, that minimizes the output power

tr[w¯¯M¯¯w¯¯](2)

with respect to the distortionless constraint

w¯¯H¯¯=I¯¯(3)

where tr[w¯¯M¯¯w¯¯] is the trace of the matrix w¯¯M¯¯w¯¯ and I¯¯ is the identity matrix. The matrix M¯¯=B¯B¯ denotes the data covariance matrix. Capon’s estimator realizing the minimal output power results in

g¯C=w¯¯B¯=[H¯¯M¯¯1H¯¯]1H¯¯M¯¯1B¯.(4)

The robustness of the method can be improved by the diagonal loading technique M¯¯M¯¯+σd2I¯¯, where σd is the so called diagonal loading parameter [3]. Thus, the estimator depends on the measurements and on the measurement positions so that measurement errors as well as measurement position errors transfer onto the estimator.

3 Error Propagation

In the following, the error propagation of Capon’s estimator is deduced. Since it is expectable that the errors are much smaller than the measurements themselves, the error propagation is derived by making use of a linearization. For the approximation of the matrix inversions that are necessary for calculating Capon’s estimator, the Neumann series bears of essential meaning for which reason we discuss it in a seperate section.

3.1 Neumann Series

The Neumann series is a special case of the functional calculus for linear operators or matrices, respectively, [14] and enables the approximation of matrix inversions.

Let T¯¯ denote a bounded matrix with norm T¯¯<1. Then,

[I¯¯T¯¯]1=k=0T¯¯k=I¯¯+T¯¯+(T¯2)(5)

where I¯¯ is the identity matrix and T¯¯k denotes the k’th power of the matrix T¯¯. The demand T¯¯<1 guarantees the convergence of the series. Thus, the Neumann series can be understood as a generalization of the geometric series for linear operators.

As a direct consequence it follows that

[I¯¯+T¯¯]1=[I¯¯(T¯¯)]1=k=0(1)kT¯¯k=I¯¯T¯¯+(T¯¯2).(6)

For the estimation of Capon’s error propagation a generalized formulation of the Neumann series is required. To expand the inverse of the sum of a matrix S¯¯ and a matrix T¯¯, where it is assumed that S¯¯ is invertible, the sum can be rewritten as

S¯¯+T¯¯=S¯¯(I¯¯+S¯¯1T¯¯).(7)

If S¯¯1T¯¯<1, the Neumann series can be applied to the sum I¯¯+S¯¯1T¯¯ resulting in

[S¯¯+T¯¯]1=[S¯¯(I¯¯+S¯¯1T¯¯)]1=(I¯¯+S¯¯1T¯¯)1S¯¯1=[k=0(1)k(S¯¯1T¯¯)k]S¯¯1=S¯¯1S¯¯1T¯¯S¯¯1+(S¯¯1T¯¯2).(8)

3.2 Measurement Errors

In the following, the error of Capon’s estimator resulting from measurement errors is derived. The model for the (temporal) statistically averaged accurate data B¯ without measurement errors is given by

B¯=H¯¯g¯+v¯,(9)

where g¯ is the wanted coefficient vector, H¯¯ denotes the shape matrix which describes the spatial dependence of the measurement positions with respect to the underlying model and v¯ denotes the (temporal) statistically averaged parts of the measurements that are not parametrized by the model H¯¯g¯ [3]. The corresponding accurate estimator resulting from Capon’s method is given by

g¯C=w¯¯B¯(10)

where

w¯¯=[H¯¯M¯¯1H¯¯]1H¯¯M¯¯1(11)

denotes the accurate filter matrix composed of the shape matrix H¯¯ and the data covariance matrix M¯¯=B¯B¯ [3].

The perturbed measurements B¯˜ can be rewritten as

B¯˜=B¯+δB¯,(12)

where δB¯ denotes the measurement error. Because of the linearity of the averaging process, the perturbed averaged measurements B¯˜ are given by

B¯˜=B¯+δB¯(13)

where δB¯ denotes the statistically averaged measurement error. The corresponding perturbed estimator results in

g˜¯C=w˜¯¯B¯˜(14)

where

w˜¯¯=[H¯¯M˜¯¯1H¯¯]1H¯¯M˜¯¯1(15)

denotes the perturbed filter matrix. Thus, the difference between the accurate and the perturbed estimator results in

g¯Cg˜¯C=w¯¯B¯w˜¯¯B¯˜=w¯¯B¯˜w¯¯δB¯w˜¯¯B¯˜=(w¯¯w˜¯¯)B¯˜w¯¯δB¯(16)

Within the practical application only the perturbed measurements B¯˜ are known and thus, only the filter matrix w˜¯¯ is known. For the calculation of the difference between the two filter matrices, in the following a linearized approximation is applied. By means of Eq. 12, the unknown accurate data covariance matrix can be rewritten as

M¯¯=B¯B¯=(B˜¯δB¯)(B˜¯δB¯)=M˜¯¯2B˜¯δB¯+δB¯δB¯(17)

where M˜¯¯=B˜¯B˜¯. We assume, that the measurement errors are much smaller than the measurements themselves, i.e., |δB¯||B˜¯|. This assumption is surely justified in the majority of applications. Considering for example the analysis of planetary magnetic fields, the measurement errors are smaller than 1nT, so that |δB¯|/|B˜¯|<1%. Thus, in the following all terms being quadratic within the errors (e.g., δB¯δB¯) will be neglected, so that the data covariance matrix results in

M¯¯=M˜¯¯2B˜¯δB¯=:M˜¯¯+ΔM¯¯.(18)

In the case of M˜¯¯1ΔM¯¯2<1, where .2 denotes the spectral norm, the application of the Neumann series (Section 3.1) delivers

M¯¯1=(M˜¯¯+ΔM¯¯)1M˜¯¯1M˜¯¯1ΔM¯¯M˜¯¯1(19)

as well as

P¯¯=[H¯¯M¯¯1H¯¯]1=(19)[H¯¯M˜¯¯1H¯¯H¯¯M˜¯¯1ΔM¯¯M˜¯¯1H¯¯]1=[H¯¯M˜¯¯1H¯¯]1+[H¯¯M˜¯¯1H¯¯]1H¯¯M˜¯¯1ΔM¯¯M˜¯¯1H¯¯[H¯¯M˜¯¯1H¯¯]1=P˜¯¯+P˜¯¯H¯¯M˜¯¯1ΔM¯¯M˜¯¯1H¯¯P˜¯¯(20)

for the unknown coefficient matrix P¯¯ [3]. Using

H¯¯M¯¯1=H¯¯M˜¯¯1H¯¯M˜¯¯1ΔM¯¯M˜¯¯1(21)

it follows that

w¯¯=P¯¯H¯¯M¯¯1=[P˜¯¯+P˜¯¯H¯¯M˜¯¯1ΔM¯¯M˜¯¯1H¯¯P˜¯¯][H¯¯M˜¯¯1H¯¯M˜¯¯1ΔM¯¯M˜¯¯1]=P˜¯¯H¯¯M˜¯¯1P˜¯¯H¯¯M˜¯¯1ΔM¯¯M˜¯¯1+P˜¯¯H¯¯M˜¯¯1ΔM¯¯M˜¯¯1H¯¯P˜¯¯H¯¯M˜¯¯1+(ΔM˜¯¯2)=w˜¯¯w˜¯¯ΔM¯¯M˜¯¯1+w˜¯¯ΔM¯¯M˜¯¯1H¯¯w˜¯¯=w˜¯¯+w˜¯¯ΔM¯¯M˜¯¯1(H¯¯w˜¯¯I¯¯),(22)

where I¯¯ again denotes the identity matrix. Thus, the difference between the filter matrices can be approximated by

w¯¯w˜¯¯=w˜¯¯ΔM¯¯M˜¯¯1(H¯¯w˜¯¯I¯¯).(23)

Inserting this approximation into Eq. 16 delivers

g¯Cg˜¯C=w˜¯¯ΔM¯¯M˜¯¯1(H¯¯w˜¯¯I¯¯)B˜¯w¯¯δB¯(24)

where

w¯¯δB¯=[w˜¯¯+w˜¯¯ΔM¯¯M˜¯¯1(H¯¯w˜¯¯I¯¯)]δB¯=w˜¯¯δB¯(25)

within the linear approximation. Further transformation for estimating the relative error of the estimator delivers

g¯Cg˜¯C=w˜¯¯ΔM¯¯M˜¯¯1(H¯¯w˜¯¯I¯¯)B˜¯w˜¯¯δB¯=w˜¯¯ΔM¯¯M˜¯¯1(H¯¯w˜¯¯B˜¯B˜¯)w˜¯¯δB¯=w˜¯¯ΔM¯¯M˜¯¯1(H¯¯g˜¯CB˜¯)w˜¯¯δB¯=w˜¯¯ΔM¯¯M˜¯¯1/2M˜¯¯1/2(H¯¯g˜¯CB˜¯)w˜¯¯δB¯.(26)

Making use of the Cauchy-Schwarz inequality yields

|g¯Cg˜¯C|2w˜¯¯ΔM¯¯M˜¯¯1/222H¯¯g˜¯CB˜¯M˜¯¯12+|w˜¯¯δB¯|2,(27)

where H¯¯g˜¯CB˜¯M˜¯¯12=|M˜¯¯1/2(H˜¯¯g˜¯CB˜¯)|2. For the calculation of Capon’s estimator, the weighted difference |M˜¯¯1/2(H˜¯¯g˜¯CB˜¯)|2 is minimized with respect to the unknown set of model parameters g˜¯, resulting in |M˜¯¯1/2(H˜¯¯g˜¯CB˜¯)|2107 [3]. A discussion about the calculation of the matrix M¯¯1/2 is given within the Appendix. Assuming that the weighted model mismatches are neglibigly small compared to the measurement errors, the deviation can be estimated upwards via

|g¯Cg˜¯C|2|w˜¯¯δB¯|2w˜¯¯22|δB¯|2,(28)

or equivalently

|g¯Cg˜¯C||g¯C|w˜¯¯2|δB¯||g¯C|=w˜¯¯2|δB¯||g¯C||B¯||B¯|(29)

for the relative error. When the measurements are adequately described by the underlying model, i.e., B¯H¯¯g¯C, it follows that

|g¯Cg˜¯C||g¯C|w˜¯¯2H¯¯2|δB¯||B¯|.(30)

Comparing this expression with the relative error of the least square fit estimator [15].

|g¯Lg˜¯L||g¯L|H¯¯+2H¯¯2|δB¯||B¯|,(31)

where H¯¯+ denotes the pseudoinverse of the shape matrix H¯¯, shows that the error propagation of Capon’s estimator follows the structure of the error propagation of the least square fit estimator. Since Capon’s method is based on the evaluation of averaged data, it can be seen that Gaussian errors with a vanishing mean value (i.e., δB¯=0) do not influence the estimator. Making use of the distortionless constraint [3].

w˜¯¯H¯¯=I¯¯(32)

delivers

w˜¯¯2=maxB¯|w˜¯¯B¯||B¯|=maxH¯¯g˜¯|w˜¯¯H¯¯g˜¯||H¯¯g˜¯|=maxg˜¯|g˜¯||H¯¯g˜¯|=(ming˜¯|H¯¯g˜¯||g˜¯|)1=1Σmin=H¯¯+2,(33)

where Σmin denotes the smallest singular value of the shape matrix H¯¯ [15], i.e. 1/Σmin describes the largest singular value of H¯¯+. Using the definition of the condition number

κ(H¯¯)=ΣmaxΣmin=H¯¯+2H¯¯2=w¯¯2H¯¯21,(34)

where Σmax denotes the largest singular value of the shape matrix H¯¯, yields

|g¯Cg˜¯C||g¯C|κ(H¯¯)|δB¯||B¯|.(35)

Thus, Capon’s estimator propagates the measurement errors in the same way as the least square fit estimator. This mathematical property can be interpreted as follows: These two methods differ in the filter matrices H¯¯+ and w¯¯ which weight or eliminate parts of the data in different ways. Since these matrices are applied to the same set of measurements, which is characterized by a given measurement error, the different weighting of the data or the elimination of subsets does not reduce the errors, since the weighted data originate from the same ensemble. The upper bound of the relative estimation error is qualitatively sketched in Figure 1.

FIGURE 1
www.frontiersin.org

FIGURE 1. Sketch of the upper bound for the relative estimation error resulting from measurement errors with respect to the condition number κ(H¯¯) of the shape matrix.

Furthermore, it should be noted that the condition number of the shape matrix determines how the measurement errors are amplified. Thereby, the condition number depends on the underlying model and the measurement positions, whereas the measurement error is produced by the sensor. For a given underlying model, the condition number solely depends on the measurement positions. Thus, the estimation errors can be reduced by analyzing measurements from suitable data points. Considering the analysis of Mercury’s magnetic field, the underlying model describes the geometry of the magnetic field, for example the internal dipole and quadrupole field [2, 3]. For the estimation of the corresponding Gauss coefficients [16, 17], the data points have to cover the geometry of the field properly. For example, the superposition of Mercury’s internal dipole field with the quadrupole field can equivalently be described as a northward shifted dipole field. When only measurement positions in the northern hemisphere are available, the condition number of the shape matrix increases, resulting in large estimation errors so that the internal field can be misinterpreted as a strong dipole field. The analysis of measurement points covering the southern and the northern hemisphere symmetrically, like that of the BepiColombo mission, decreases the condition number and enables a more detailed characterization of the geometry of Mercury’s internal magnetic field [18]. The smaller the condition number of the shape matrix becomes, the better the geometry of the field is covered by the data points [18].

Before discussing the errors resulting from the perturbed determination of the measurement positions, let us make a general comment about the error of Capon’s estimator resulting from measurement errors. The perturbed data can be rewritten in the form

B˜¯=B¯+δB¯=H¯¯g¯+v¯+δB¯=H¯¯g¯+v˜¯,(36)

where v˜¯=v¯+δB¯. Since the filter matrix w˜¯¯ eliminates the non-parametrized parts from the total measured field [3], one might suppose that the measurement errors do not influence the estimation since these errors can as well be interpreted as non-parametrized parts. From the elimination of the non-parametrized parts

0=w˜¯¯v˜¯=w˜¯¯v¯+w˜¯¯δB¯(37)

it does not follow that w˜¯¯δB¯=0. Since this term dominates the error of the estimator (cf. Eq. 28), the measurement errors are not eliminated by the filter matrix in general.

3.3 Measurement Position Errors

The model for the correctly determined measurement positions is given by

B¯=H¯¯g¯+v¯.(38)

The corresponding accurate estimator results in

g¯C=w¯¯B¯,(39)

where

w¯¯=[H¯¯M¯¯1H¯¯]1H¯¯M¯¯1.(40)

The perturbed determination of the sensor’s positions transfers onto a perturbed shape matrix H¯¯. Thus, the underlying model H¯¯g¯ and the non-parametrized parts v¯ are perturbed so that the noisy model is given by

B¯=H˜¯¯g˜¯+v˜¯(41)

and the corresponding perturbed estimator results in

g˜¯C=w˜¯¯B¯,(42)

where

w˜¯¯=[H˜¯¯M¯¯1H˜¯¯]1H˜¯¯M¯¯1.(43)

The deviation between the accurate and the perturbed estimator results in

g¯Cg˜¯C=(w¯¯w˜¯¯)B¯.(44)

As discussed above, only the perturbed filter matrix H˜¯¯ is known. Within a linear approximation, the unknown matrix H¯¯ can be rewritten as H¯¯=H˜¯¯+ΔH¯¯, where ΔH¯¯2H˜¯¯2. For example, consider the shape matrix

H¯¯=[1(RMr1)31(RMr2)3](45)

which describes the magnetic field of a current sheet superposed with Mercury’s internal dipolar field [6], where RM indicates the planetary radius of Mercury and r1 and r2 are the measurement positions. The perturbed filter matrix is given by

H˜¯¯=[1(RMr1+Δr1)31(RMr2+Δr2)3]=[1RM3r13(1+Δr1/r1)31RM3r23(1+Δr2/r2)3].(46)

Assuming that Δriri, for i=1,2, the perturbed matrix can be rewritten as

H˜¯¯=[1(RMr1)31(RMr2)3]+[03RM3r13Δr1r103RM3r23Δr2r2]=H¯¯ΔH¯¯(47)

by making use of a Taylor series expansion. It should be noted that the shape matrix can be linearized for any underlying model by performing a Taylor series expansion for the functions within the model. Linearization of the error terms for the estimation of the unknown accurate filter matrix delivers

w¯¯=(40)[H¯¯M¯¯1H¯¯]1H¯¯M¯¯1=[(H˜¯¯+ΔH¯¯)M¯¯1(H˜¯¯+ΔH¯¯)]1(H˜¯¯+ΔH¯¯)M¯¯1=[H˜¯¯M¯¯1(H˜¯¯+ΔH¯¯)+ΔH¯¯M¯¯1(H˜¯¯+ΔH¯¯)]1(H˜¯¯+ΔH¯¯)M¯¯1=[H˜¯¯M¯¯1H˜¯¯+H˜¯¯M¯¯1ΔH¯¯+ΔH¯¯M¯¯1H˜¯¯]1(H˜¯¯+ΔH¯¯)M¯¯1.(48)

Making use of the Neumann series results in

[H˜¯¯M¯¯1H˜¯¯+H˜¯¯M¯¯1ΔH¯¯+ΔH¯¯M¯¯1H˜¯¯]1=[H˜¯¯M¯¯1H˜¯¯]1[H˜¯¯M¯¯1H˜¯¯]1[H˜¯¯M¯¯1ΔH¯¯+ΔH¯¯M¯¯1H˜¯¯][H˜¯¯M¯¯1H˜¯¯]1=P˜¯¯P˜¯¯H˜¯¯M¯¯1ΔH¯¯P˜¯¯P˜¯¯ΔH¯¯M¯¯1H˜¯¯P˜¯¯,(49)

so that

w¯¯=P˜¯¯H˜¯¯M¯¯1P˜¯¯H˜¯¯M¯¯1ΔH¯¯P˜¯¯H˜¯¯M¯¯1P˜¯¯ΔH¯¯M¯¯1H˜¯¯P˜¯¯H˜¯¯M¯¯1+P˜¯¯ΔH¯¯M¯¯1=w˜¯¯w˜¯¯ΔH¯¯w˜¯¯P˜¯¯ΔH¯¯M¯¯1H˜¯¯w˜¯¯+P˜¯¯ΔH¯¯M¯¯1=w˜¯¯w˜¯¯ΔH¯¯w˜¯¯P˜¯¯ΔH¯¯M¯¯1(H˜¯¯w˜¯¯I¯¯),(50)

or equivalently

w¯¯w˜¯¯=w˜¯¯ΔH¯¯w˜¯¯P˜¯¯ΔH¯¯M¯¯1(H˜¯¯w˜¯¯I¯¯).(51)

Thus, the deviation of the estimator can be approximated via

g¯Cg˜¯C=w˜¯¯ΔH¯¯w˜¯¯B¯P˜¯¯ΔH¯¯M¯¯1(H˜¯¯w˜¯¯B¯B¯).(52)

By w˜¯¯B¯=g˜¯C it follows that

g¯Cg˜¯C=w˜¯¯ΔH¯¯g˜¯CP˜¯¯ΔH¯¯M¯¯1(H˜¯¯g˜¯CB¯),(53)

and therefore

|g¯Cg˜¯C|2w˜¯¯22ΔH¯¯22|g˜¯C|2+P˜¯¯ΔH¯¯M¯¯1/222|M¯¯1/2(H˜¯¯g˜¯CB¯)|2.(54)

Within the application of Capon’s method to Mercury’s magnetic field analysis, the first summand on the right hand side of Eq. 54 is of the order of about 10nT2, whereas P˜¯¯ΔH¯¯M¯¯1/222|M¯¯1/2(H˜¯¯g˜¯CB¯)|2500107nT2 and thus, the weighted model mismatches are negligibly small compared to the measurement errors, so that the deviation can be estimated upwards via

|g¯Cg˜¯C|2w˜¯¯22ΔH¯¯22|g˜¯C|2.(55)

The relative error results in

|g¯Cg˜¯C||g˜¯C|w˜¯¯2H˜¯¯2ΔH¯¯2H˜¯¯2=κ(H˜¯¯)ΔH¯¯2H˜¯¯2,(56)

where κ(H˜¯¯) denotes the condition number of the perturbed shape matrix H˜¯¯. The upper bound for the relative error with respect to the condition number is qualitatively sketched in Figure 2. Since w˜¯¯2=H˜¯¯+2 (cf. Eq. 33), the error propagation of Capon’s estimator again equals the error propagation of the least square fit estimator [15].

|g¯Lg˜¯L||g˜¯L|H˜¯¯+2H˜¯¯2ΔH¯¯2H˜¯¯2=κ(H˜¯¯)ΔH¯¯2H˜¯¯2.(57)

FIGURE 2
www.frontiersin.org

FIGURE 2. Sketch of the upper bound for the relative estimation error resulting from measurement position errors with respect to the condition number κ(H˜¯¯) of the shape matrix.

3.4 Measurement Errors and Measurement Position Errors

Within the former sections the influences of measurement errors and measurement position errors have been discussed separately. Within the practical application it is expectable that both errors occur simultanously. The relative error resulting from the noisy measurements and the perturbed measurement positions is given by the quadratic sum of the two cases discussed above

|g¯Cg˜¯C|2|g˜¯C|2κ2(H˜¯¯)(ΔH¯¯22H˜¯¯22+|δB¯|2|B˜¯|2).(58)

This can be derived as follows:

The accurate model is given by

B¯=H¯¯g¯+v¯,(59)

so that the corresponding accurate estimator results in

g¯C=w¯¯B¯(60)

where

w¯¯=[H¯¯M¯¯1H¯¯]1H¯¯M¯¯1.(61)

The measured data are not affected by the perturbed determination of the sensor’s positions. The data are solely affected by measurement errors δB¯ so that the perturbed model can be written as

B˜¯=B¯+δB¯=H˜¯¯g˜¯+v˜¯+δB¯(62)

The corresponding perturbed estimator results in

g˜¯C=w˜¯¯B˜¯,(63)

where

w˜¯¯=[H˜¯¯M˜¯¯1H˜¯¯]1H˜¯¯M˜¯¯1(64)

and M¯¯=M˜¯¯+ΔM¯¯. The noisy shape matrix can again be written as H¯¯=H˜¯¯+ΔH¯¯ within a linear approximation. Thus, the deviation between the accurate and the perturbed estimator is given by

g¯Cg˜¯C=w¯¯B¯w˜¯¯B˜¯=w¯¯B˜¯w¯¯δB¯w˜¯¯B˜¯=(w¯¯w˜¯¯)B˜¯w¯¯δB¯.(65)

By making use of the Neumann series, the unknown filter matrix w¯¯ can be rewritten within a linear approximation resulting in

w¯¯=[H¯¯M¯¯1H¯¯]1H¯¯M¯¯1=[(H˜¯¯+ΔH¯¯)(M˜¯¯+ΔM¯¯)1(H˜¯¯+ΔH¯¯)]1(H˜¯¯+ΔH¯¯)(M˜¯¯+ΔM¯¯)1=[(H˜¯¯+ΔH¯¯)(M˜¯¯1M˜¯¯1ΔM¯¯M˜¯¯1)(H˜¯¯+ΔH¯¯)]1(H˜¯¯+ΔH¯¯)(M˜¯¯1M˜¯¯1ΔM¯¯M˜¯¯1).

Using

(H˜¯¯+ΔH¯¯)(M˜¯¯1M˜¯¯1ΔM¯¯M˜¯¯1)(H˜¯¯+ΔH¯¯)=(H˜¯¯M˜¯¯1H˜¯¯M˜¯¯1ΔM¯¯M˜¯¯1+ΔH¯¯M˜¯¯1)(H˜¯¯+ΔH¯¯)=H˜¯¯M˜¯¯1H˜¯¯H˜¯¯M˜¯¯1ΔM¯¯M˜¯¯1H˜¯¯+ΔH¯¯M˜¯¯1H˜¯¯+H˜¯¯M˜¯¯1ΔH¯¯=P˜¯¯1H˜¯¯M˜¯¯1ΔM¯¯M˜¯¯1H˜¯¯+ΔH¯¯M˜¯¯1H˜¯¯+H˜¯¯M˜¯¯1ΔH¯¯(66)

and again making use of the Neumann series delivers

w¯¯=[P˜¯¯+P˜¯¯H˜¯¯M˜¯¯1ΔM¯¯M˜¯¯1H˜¯¯P˜¯¯P˜¯¯ΔH¯¯M˜¯¯1H˜¯¯P˜¯¯P˜¯¯H˜¯¯M˜¯¯1ΔH¯¯P˜¯¯](H˜¯¯M˜¯¯1H˜¯¯M˜¯¯1ΔM¯¯M˜¯¯1+ΔH¯¯M˜¯¯1)=P˜¯¯H˜¯¯M˜¯¯1P˜¯¯H˜¯¯M˜¯¯1ΔM¯¯M˜¯¯1+P˜¯¯ΔH¯¯M˜¯¯1+P˜¯¯H˜¯¯M˜¯¯1ΔM¯¯M˜¯¯1H˜¯¯P˜¯¯H˜¯¯M˜¯¯1P˜¯¯ΔH¯¯M˜¯¯1H˜¯¯P˜¯¯H˜¯¯M˜¯¯1P˜¯¯H˜¯¯M˜¯¯1ΔH¯¯P˜¯¯H˜¯¯M˜¯¯1=w˜¯¯w˜¯¯ΔM¯¯M˜¯¯1+P˜¯¯ΔH¯¯M˜¯¯1+w˜¯¯ΔM¯¯M˜¯¯1H˜¯¯w˜¯¯P˜¯¯ΔH¯¯M˜¯¯1H˜¯¯w˜¯¯w˜¯¯ΔH¯¯w˜¯¯.(67)

Taking into account that

w¯¯δB¯=w˜¯¯δB¯(68)

within the linear approximation as well as

w¯¯w˜¯¯=w˜¯¯ΔM¯¯M˜¯¯1+P˜¯¯ΔH¯¯M˜¯¯1+w˜¯¯ΔM¯¯M˜¯¯1H˜¯¯w˜¯¯P˜¯¯ΔH¯¯M˜¯¯1H˜¯¯w˜¯¯w˜¯¯ΔH¯¯w˜¯¯,(69)

the deviation between the estimators can be approximated by

g¯Cg˜¯C=(w¯¯w˜¯¯)B˜¯w¯¯δB¯=w˜¯¯ΔM¯¯M˜¯¯1B˜¯+P˜¯¯ΔH¯¯M˜¯¯1B˜¯+w˜¯¯ΔM¯¯M˜¯¯1H˜¯¯w˜¯¯B˜¯P˜¯¯ΔH¯¯M˜¯¯1H˜¯¯w˜¯¯B˜¯w˜¯¯ΔH¯¯w˜¯¯B˜¯w˜¯¯δB¯(70)

Using w˜¯¯B˜¯=g˜¯C, it follows that

g¯Cg˜¯C=w˜¯¯ΔM¯¯M˜¯¯1B˜¯+P˜¯¯ΔH¯¯M˜¯¯1B˜¯+w˜¯¯ΔM¯¯M˜¯¯1H˜¯¯g˜¯CP˜¯¯ΔH¯¯M˜¯¯1H˜¯¯g˜¯Cw˜¯¯ΔH¯¯g˜¯Cw˜¯¯δB¯=w˜¯¯ΔM¯¯M˜¯¯1(H˜¯¯g˜¯CB˜¯)+P˜¯¯ΔH¯¯M˜¯¯1(B˜¯H˜¯¯g˜¯C)w˜¯¯ΔH¯¯g˜¯Cw˜¯¯δB¯(71)

and thus,

|g¯Cg˜¯C|2=w˜¯¯ΔM¯¯M˜¯¯1/2M˜¯¯1/2(H˜¯¯g˜¯CB˜¯)+P˜¯¯ΔH¯¯M˜¯¯1/2M˜¯¯1/2(B˜¯H˜¯¯g˜¯C)w˜¯¯ΔH¯¯g˜¯Cw˜¯¯δB¯22w˜¯¯ΔM¯¯M˜¯¯1/222|M˜¯¯1/2(H˜¯¯g˜¯CB˜¯)|2+P˜¯¯ΔH¯¯M˜¯¯1/222|M˜¯¯1/2(B˜¯H˜¯¯g˜¯C)|2+w˜¯¯22ΔH¯¯22|g˜¯C|2+w˜¯¯22|δB¯|2w˜¯¯22ΔH¯¯22|g˜¯C|2+w˜¯¯22|δB¯|2.(72)

Assuming that B˜¯H˜¯¯g˜¯C the relative error can be estimated by

|g¯Cg˜¯C|2|g˜¯C|2=w˜¯¯22ΔH¯¯22+w˜¯¯22|δB¯|2|B˜¯|2|B˜¯|2|g˜¯C|2=w˜¯¯22H˜¯¯22ΔH¯¯22H˜¯¯22+w˜¯¯22|δB¯|2|B˜¯|2|H˜¯¯g˜¯C|2|g˜¯C|2w˜¯¯22H˜¯¯22ΔH¯¯22H˜¯¯22+w˜¯¯22||H˜¯¯||22|δB¯|2|B˜¯|2=w˜¯¯22H˜¯¯22(ΔH¯¯22H˜¯¯22+|δB¯|2|B˜¯|2)=κ2(H˜¯¯)(ΔH¯¯22H˜¯¯22+|δB¯|2|B˜¯|2)(73)

It is important to note that the right sides of Eqs. 35, 56 and 73 represent upper bounds for the error of the estimator. These bounds can be much larger than the true errors. The great advantage of the bounds lies in the fact that they solely depend on known quantities, wheras the accurate estimator for calculating the true estimation error is unknown within the practical application of the method. Thus, the above derived upper bounds are conservative and guarantee that the true estimation errors cannot exceed these bounds as long as the measurement errors as well as the measurement position errors are much smaller than the measurements themselves so that the linear approximation is valid.

4 Finite Sample Averaging

For the calculation of Capon’s estimator, averaged magnetic field data

B¯=1Qα=1QB¯α(74)

are required [3]. Here, Q denotes the number of measurements at a fixed set of data points. Within the practical application of the method only a finite number Q< of samples is available, which yields a standard deviation of Capon’s estimator g¯C. In the following, an approximation for the variance of Capon’s estimator is derived by regarding Capon’s method as a special case of the maximum likelihood method [6].

In the vicinity of its maximum value the likelihood function can be approximated as

exp[12α=1Q(B¯αH¯¯g¯)N¯¯1(B¯αH¯¯g¯)],(75)

where N¯¯ denotes the noise covariance matrix [6, 19]. Before deducing the error resulting from the finite number of samples let us make a general comment about the construction of the likelihood function in terms of the averaging process. The noise matrix is already statistically evaluated so that the likelihood function is constructed by making use of two different averaging processes. Within the first averaging process the noise matrix N¯¯ is calculated. This averaging process does not incorporate the underlying model since only the distribution of the data around the mean value B¯ is determined which allows to assess the general quality of the measurements. When the distribution of the measurements is large, significant errors within the subsequent data fitting are expectable. Especially, these errors are independent of the chosen model. Within the second averaging process a specific underlying model H¯¯g¯ is fitted to the data B¯α.

For most practical applications, the noise matrix N¯¯ is unknown and has to be approximated. In the special case of Gaussian errors with zero mean and variance σn, the noise covariance matrix can be written as N¯¯=σn2I¯¯, where I¯¯ denotes the identity matrix. Substituting the noise covariance matrix N¯¯ by the data covariance matrix M¯¯=B¯B¯, the maximum likelihood estimator may be converted into Capon’s estimator [6]. The corresponding likelihood function is modified to

exp[12α=1Q(B¯αH¯¯g¯)M¯¯1(B¯αH¯¯g¯)].(76)

The weighted difference between the model and the data can be rewritten as

(B¯αH¯¯g¯)M¯¯1(B¯αH¯¯g¯)=(B¯αg¯H¯¯)M¯¯1(B¯αH¯¯g¯)=B¯αM¯¯1B¯α2g¯H¯¯M¯¯1B¯α+g¯H¯¯M¯¯1H¯¯g¯.(77)

As discussed above, the data covariance matrix M¯¯ is already statistically evaluated so that the averaging process results in

α=1Q(B¯αH¯¯g¯)M¯¯1(B¯αH¯¯g¯)=QB¯M¯¯1B¯2Qg¯H¯¯M¯¯1B¯+Qg¯H¯¯M¯¯1H¯¯g¯.(78)

Insertion into Eq. 76 yields

exp{Q2[B¯M¯¯1B¯2g¯H¯¯M¯¯1B¯+g¯H¯¯M¯¯1H¯¯g¯]}(79)

or equivalently

exp{Q2[BiMij1Bj2gkHkiMij1Bj+glHlkMki1Hijgj]}.(80)

The m’th component of Capon’s estimator g¯C corresponds with the maximum value of the likelihood function [6, 19]. The corresponding standard deviation (1σ-error) is described by the width of the likelihood function at its maximum value which is given by [6, 19].

σgm=[gm2ln]1/2,(81)

where

gm2ln=gm(1gm).(82)

Since

gm=Q2gm[BiMij1Bj2gkHkiMij1Bj+glHlkMki1Hijgj]=Q2[2δkmHkiMij1Bj+δlmHlkMki1Hijgj+glHlkMki1Hijδjm]=Q2[2HmiMij1Bj+HmkMki1Hijgj+glHlkMki1Him](83)

where δij denotes the Kronecker delta, the m’th component of Capon’s estimator g¯C can be calculated via [6].

gm=0(84)

yielding

2HmiMij1Bj+HmkMki1Hijgj+glHlkMki1Him=0.(85)

Because of Mki1=Mik1, as well as Hlk=Hkl and Him=Hmi in the case of real shape matrices [3], it follows that

0=2HmiMij1Bj+HmkMki1Hijgj+glHklMik1Hmi=2HmiMij1Bj+HmkMki1Hijgj+HmiMik1Hklgl=2HmiMij1Bj+2HmkMki1Hijgj(86)

and therefore,

g¯C=[H¯¯M¯¯1H¯¯]1H¯¯M¯¯1B¯(87)

which is in agreement with the estimator resulting from the linear algebraic formulation of Capon’s method [3]. The second derivative results in

gm2ln=Q2gm[2HmiMij1Bj+HmkMki1Hijgj+glHlkMki1Him]=Q2[HmkMki1Hijδjm+δlmHlkMki1Him]=Q2[HmkMki1Him+HmkMki1Him]=QHmkMki1Him.(88)

Using the definition of the coefficient matrix [3].

P¯¯=[H¯¯M¯¯1H¯¯]1(89)

or equivalently

(P¯¯1)mm=HmkMki1Him(90)

delivers

gm2ln=Q(P¯¯1)mm.(91)

Thus, the error of the m’th component of Capon’s estimator results in

σgm=1Q(P¯¯1)mm,(92)

which declines as 1/Q. The functional dependence of the error is qualitatively illustrated in Figure 3.

FIGURE 3
www.frontiersin.org

FIGURE 3. Sketch of the m’the component of Capon’s estimator and the corresponding 1σ-error subject to the sample size Q. The error declines as 1/Q.

5 Summary and Outlook

The analysis of the error propagation is of major importance for the application of linear inversion methods. Within a linear approximation, upper bounds for the errors of Capon’s estimator resulting from measurement errors and measurement position errors are derived. These upper bounds solely depend on known quantities, i.e., measurements and measurement positions, whereas the true estimation error cannot be calculated within the practical application of the method, since the accurate estimator is unavailable. It turns out that Capon’s method provides the same error propagation as the least square fit method. These two methods differ in the filter matrices which weight or eliminate parts of the data in different ways. Since these matrices are applied to the same set of measurements, the different weighting of the data or the elimination of subsets does not reduce the errors. The condition number of the shape matrix is the key parameter for the error propagation, since it determines how the errors are amplified. The measurement errors as well as the measurement position errors have to be estimated from the measurements. For a given underlying model, the condition number of the shape matrix solely depends on the measurement positions. Thus, the amplification of the errors can be reduced by choosing preferred data points.

Furthermore, Capon’s method is based on the evaluation of statistically averaged data. For the practical application of the method only a finite number Q of samples is available. This limited number of samples results in an error of Capon’s estimator which can be derived by regarding Capon’s method as a special case of the maximum likelihood estimator. The more samples are available, the smaller the error becomes, since the error declines as 1/Q.

As a follow-up of the generalized derivation of Capon’s method, the present work establishes the mathematical basis of Capon’s error propagation for the practical application of the method.

Appendix: Matrix Power of the Data Covariance Matrix

For the calculation of Capon’s estimator the inverse data covariance matrix M¯¯1 is necessary. First of all it should be noted that the inverse exists due to the averaging process, whereas the matrix B¯B¯ is not invertible [3]. In the case of a vanishing standard deviation (σn=0), the application of the diagonal loading technique [3] guarantees the existence of the inverse data covariance matrix, since the condition number of M¯¯ is close to unity at the optimal diagonal loading parameter.

Furthermore, there exists a unitary transformation V¯¯ so that

M¯¯=V¯¯D¯¯M¯¯V¯¯1(93)

where

D¯¯M¯¯=diag[λ1,,λn](94)

is a diagonal matrix that contains the eigenvalues of M¯¯. Therefore, the inverse data covariance matrix results in

M¯¯1=[V¯¯D¯¯M¯¯V¯¯1]1=V¯¯D¯¯M¯¯1V¯¯1,(95)

where

D¯¯M¯¯1=diag[λ11,,λn1].(96)

Further transformation delivers

M¯¯1=V¯¯D¯¯M¯¯1V¯¯1=V¯¯D¯¯M¯¯1/2D¯¯M¯¯1/2V¯¯1=V¯¯D¯¯M¯¯1/2V¯¯1V¯¯D¯¯M¯¯1/2V¯¯1=M¯¯1/2M¯¯1/2(97)

so that the square root of the inverse data covariance matrix is defined as

M¯¯1/2=V¯¯D¯¯M¯¯1/2V¯¯1(98)

where

D¯¯M¯¯1/2=diag[λ11/2,,λn1/2].(99)

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

All authors contributed conception and design of the study; ST and UM wrote the first draft of the manuscript; All authors contributed to manuscript revision, read and approved the submitted version.

Funding

We acknowledge support by the German Research Foundation and the Open Access Publication Funds of the Technische Universität Braunschweig. The work by YN is supported by the Austrian Space Applications Program at the Austrian Research Promotion Agency under contract 865967. DH was supported by the German Ministerium für Wirtschaft und Energie and the German Zentrum für Luft-und Raumfahrt under contract 50 QW1501.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer LZ declared a past collaboration with one of the authors YN to the handling editor.

Acknowledgments

The authors are grateful for stimulating discussions and helpful suggestions by Karl-Heinz Glassmeier, Patrick Kolhey and Alexander Schwenke.

References

1. Capon J. High-resolution Frequency-Wavenumber Spectrum Analysis. Proc IEEE (1969) 57:1408–18. doi:10.1109/PROC.1969.7278

CrossRef Full Text | Google Scholar

2. Toepfer S, Narita Y, Heyner D, Motschmann U. The Capon Method for Mercury's Magnetic Field Analysis. Front Phys (2020) 8:249. doi:10.3389/fphy.2020.00249

CrossRef Full Text | Google Scholar

3. Toepfer S, Narita Y, Heyner D, Kolhey P, Motschmann U. Mathematical Foundation of Capon's Method for Planetary Magnetic Field Analysis. Geosci Instrum Method Data Syst (2020) 9:471–81. doi:10.5194/gi-9-471-2020

CrossRef Full Text | Google Scholar

4. Motschmann U, Woodward TI, Glassmeier KH, Southwood DJ, Pinçon JL. Wavelength and Direction Filtering by Magnetic Measurements at Satellite Arrays: Generalized Minimum Variance Analysis. J Geophys Res (1996) 101:4961–5. doi:10.1029/95JA03471

CrossRef Full Text | Google Scholar

5. Su W, Gu H, Ni J, Liu G. Performance Analysis of MVDR Algorithm in the Presence of Amplitude and Phase Errors. IEEE Trans Antennas Propagat (2001) 49(12):1875–7. doi:10.1109/8.982472

CrossRef Full Text | Google Scholar

6. Narita Y. A Note on Capon’s Minimum Variance Projection for Multi-Spacecraft Data Analysis. Front Phys (2019) 7:121. doi:10.3398/fphy.2019.00008

CrossRef Full Text | Google Scholar

7. Sahraoui F, Belmont G, Rezeau L, Cornilleau-Wehrlin N, Pinçon JL, Balogh A. Anisotropic Turbulent Spectra in the Terrestrial Magnetosheath as Seen by the Cluster Spacecraft. Phys Rev Lett (2006) 96:7. doi:10.1103/PhysRevLett.96.075002

CrossRef Full Text | Google Scholar

8. Huang SY, Zhou M, Sahraoui F, Deng XH, Pang Y, Yuan ZG, et al. Wave Properties in the Magnetic Reconnection Diffusion Region with Highβ: Application of Thek-Filtering Method to Cluster Multispacecraft Data. J Geophys Res (2010) 115:17–9. doi:10.1029/2010JA015335

CrossRef Full Text | Google Scholar

9. Narita Y, Comişel H, Motschmann U. Spatial Structure of Ion-Scale Plasma Turbulence. Front Phys (2014) 2:13. doi:10.3389/fphy.2014.00013

CrossRef Full Text | Google Scholar

10. Roberts OW, Li X, Jeska L. A Statistical Study of the Solar Wind Turbulence at Ion Kinetic Scales Using the K-Filtering Technique and Cluster Data. ApJ (2015) 802:1. doi:10.1088/0004-637X/802/1/2

CrossRef Full Text | Google Scholar

11. Wang T, Cao J, Fu H, Meng X, Dunlop M. Compressible Turbulence with Slow-Mode Waves Observed in the Bursty Bulk Flow of Plasma Sheet. Geophys Res Lett (2016) 43(5):1854–61. doi:10.1002/2016GL068147

CrossRef Full Text | Google Scholar

12. Zhang L, He J, Narita Y, Feng X. Reconstruction Test of Turbulence Power Spectra in 3D Wavenumber Space with at Most 9 Virtual Spacecraft Measurements. J Geophys Res Space Phys (2021) 126:e2019JA027413. doi:10.1029/2019JA027413

CrossRef Full Text | Google Scholar

13. Narita Y, Plaschke F, Magnes W, Fischer D, Schmid D. Error Estimate for Fluxgate Magnetometer In-Flight Calibration on a Spinning Spacecraft. Geosci Instrum Method Data Syst (2021) 10:13–24. doi:10.5194/gi-10-13-202

CrossRef Full Text | Google Scholar

14. Werner D. Funktionalanalysis. Berlin: Springer (2007)

15. Belsley DA, Kuh E, Welsch RE. The Condition Number, Regression Diagnostics: Identifying Influential Data and Sources of Collinearity. New York: Wiley (1980). p. 100–4.

16. Gauß CF. (1839) Allgemeine Theorie des Erdmagnetismus: Resultate aus den Beobachtungen des magnetischen Vereins im Jahre 1838, editorsCF Gauss, and W Weber. Leipzig: Weidmannsche Buchhandlung (1839). p. 1–57.

17. Glassmeier K-H, Tsurutani BT. Carl Friedrich Gauss—General Theory of Terrestrial Magnetism—A Revised Translation of the German Text. Hist Geo Space Sci (2014) 5:11–62. doi:10.5194/hgss-5-11-2014

CrossRef Full Text | Google Scholar

18. Toepfer S, Narita Y, Glassmeier K-H, Heyner D, Kolhey P, Motschmann U, et al. The Mie Representation for Mercury's Magnetic Field. Earth Planets Space (2021) 73:65. doi:10.1186/s40623-021-01386-4

CrossRef Full Text | Google Scholar

19. Dodelson S. Modern Cosmology. London: Academic Press (2003).

Keywords: Capon’s method, error propagation, least-squares method, maximum likelihood, condition number

Citation: Toepfer S, Narita Y, Heyner D and Motschmann U (2021) Error Propagation of Capon’s Minimum Variance Estimator. Front. Phys. 9:684410. doi: 10.3389/fphy.2021.684410

Received: 23 March 2021; Accepted: 20 May 2021;
Published: 09 June 2021.

Edited by:

Jiansen He, Peking University, China

Reviewed by:

Tieyan Wang, Rutherford Appleton Laboratory, United Kingdom
Lei Zhang, China Academy of Space Technology (CAST), China

Copyright © 2021 Toepfer, Narita, Heyner and Motschmann. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: S. Toepfer, s.toepfer@tu-braunschweig.de

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.