Abstract
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 [1–3], 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 [7–12]. 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 into parametrized parts , where contains the corresponding wanted model coefficients and the shape matrix describes the distribution of the measurement positions with respect to the underlying model, non-parametrized parts as well as measurement noise , so thatis valid. The measurement noise is assumed to be Gaussian with variance and zero mean . Since the shape matrix is not invertible and the non-parametrized parts are unknown, the exact solution for the wanted model coefficients is not available in general. Capon’s method delivers an estimator for the ideal solution . The method is based on the construction of a filter matrix , that minimizes the output powerwith respect to the distortionless constraintwhere is the trace of the matrix and is the identity matrix. The matrix denotes the data covariance matrix. Capon’s estimator realizing the minimal output power results inThe robustness of the method can be improved by the diagonal loading technique , where 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 denote a bounded matrix with norm . Then,where is the identity matrix and denotes the k’th power of the matrix . The demand 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 thatFor 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 and a matrix , where it is assumed that is invertible, the sum can be rewritten asIf , the Neumann series can be applied to the sum resulting in
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 without measurement errors is given bywhere is the wanted coefficient vector, denotes the shape matrix which describes the spatial dependence of the measurement positions with respect to the underlying model and denotes the (temporal) statistically averaged parts of the measurements that are not parametrized by the model [3]. The corresponding accurate estimator resulting from Capon’s method is given bywheredenotes the accurate filter matrix composed of the shape matrix and the data covariance matrix [3].
The perturbed measurements can be rewritten aswhere denotes the measurement error. Because of the linearity of the averaging process, the perturbed averaged measurements are given bywhere denotes the statistically averaged measurement error. The corresponding perturbed estimator results inwheredenotes the perturbed filter matrix. Thus, the difference between the accurate and the perturbed estimator results inWithin the practical application only the perturbed measurements are known and thus, only the filter matrix 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 aswhere . We assume, that the measurement errors are much smaller than the measurements themselves, i.e., . 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 , so that . Thus, in the following all terms being quadratic within the errors (e.g., ) will be neglected, so that the data covariance matrix results inIn the case of , where denotes the spectral norm, the application of the Neumann series (Section 3.1) deliversas well asfor the unknown coefficient matrix [3]. Usingit follows thatwhere again denotes the identity matrix. Thus, the difference between the filter matrices can be approximated byInserting this approximation into Eq. 16 deliverswherewithin the linear approximation. Further transformation for estimating the relative error of the estimator deliversMaking use of the Cauchy-Schwarz inequality yieldswhere . For the calculation of Capon’s estimator, the weighted difference is minimized with respect to the unknown set of model parameters , resulting in [3]. A discussion about the calculation of the matrix 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 viaor equivalentlyfor the relative error. When the measurements are adequately described by the underlying model, i.e., , it follows thatComparing this expression with the relative error of the least square fit estimator [15].where denotes the pseudoinverse of the shape matrix , 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., ) do not influence the estimator. Making use of the distortionless constraint [3].deliverswhere denotes the smallest singular value of the shape matrix [15], i.e. describes the largest singular value of . Using the definition of the condition numberwhere denotes the largest singular value of the shape matrix , yieldsThus, 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 and 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
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 formwhere . Since the filter matrix 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 partsit does not follow that . 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 byThe corresponding accurate estimator results inwhereThe perturbed determination of the sensor’s positions transfers onto a perturbed shape matrix . Thus, the underlying model and the non-parametrized parts are perturbed so that the noisy model is given byand the corresponding perturbed estimator results inwhereThe deviation between the accurate and the perturbed estimator results inAs discussed above, only the perturbed filter matrix is known. Within a linear approximation, the unknown matrix can be rewritten as , where . For example, consider the shape matrixwhich describes the magnetic field of a current sheet superposed with Mercury’s internal dipolar field [6], where indicates the planetary radius of Mercury and and are the measurement positions. The perturbed filter matrix is given byAssuming that , for , the perturbed matrix can be rewritten asby 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 deliversMaking use of the Neumann series results inso thator equivalentlyThus, the deviation of the estimator can be approximated viaBy it follows thatand thereforeWithin 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 , whereas and thus, the weighted model mismatches are negligibly small compared to the measurement errors, so that the deviation can be estimated upwards viaThe relative error results inwhere denotes the condition number of the perturbed shape matrix . The upper bound for the relative error with respect to the condition number is qualitatively sketched in Figure 2. Since (cf. Eq. 33), the error propagation of Capon’s estimator again equals the error propagation of the least square fit estimator [15].
FIGURE 2
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 aboveThis can be derived as follows:
The accurate model is given byso that the corresponding accurate estimator results inwhereThe measured data are not affected by the perturbed determination of the sensor’s positions. The data are solely affected by measurement errors so that the perturbed model can be written asThe corresponding perturbed estimator results inwhereand . The noisy shape matrix can again be written as within a linear approximation. Thus, the deviation between the accurate and the perturbed estimator is given byBy making use of the Neumann series, the unknown filter matrix can be rewritten within a linear approximation resulting inUsingand again making use of the Neumann series deliversTaking into account thatwithin the linear approximation as well asthe deviation between the estimators can be approximated byUsing , it follows thatand thus,Assuming that the relative error can be estimated byIt 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 dataare 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 of samples is available, which yields a standard deviation of Capon’s estimator . 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 aswhere 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 is calculated. This averaging process does not incorporate the underlying model since only the distribution of the data around the mean value 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 is fitted to the data .
For most practical applications, the noise matrix is unknown and has to be approximated. In the special case of Gaussian errors with zero mean and variance , the noise covariance matrix can be written as , where denotes the identity matrix. Substituting the noise covariance matrix by the data covariance matrix , the maximum likelihood estimator may be converted into Capon’s estimator [6]. The corresponding likelihood function is modified toThe weighted difference between the model and the data can be rewritten asAs discussed above, the data covariance matrix is already statistically evaluated so that the averaging process results inInsertion into Eq. 76 yieldsor equivalentlyThe m’th component of Capon’s estimator corresponds with the maximum value of the likelihood function [6, 19]. The corresponding standard deviation (-error) is described by the width of the likelihood function at its maximum value which is given by [6, 19].whereSincewhere denotes the Kronecker delta, the m’th component of Capon’s estimator can be calculated via [6].yieldingBecause of , as well as and in the case of real shape matrices [3], it follows thatand therefore,which is in agreement with the estimator resulting from the linear algebraic formulation of Capon’s method [3]. The second derivative results inUsing the definition of the coefficient matrix [3].or equivalentlydeliversThus, the error of the m’th component of Capon’s estimator results inwhich declines as . The functional dependence of the error is qualitatively illustrated in Figure 3.
FIGURE 3
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 .
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 is necessary. First of all it should be noted that the inverse exists due to the averaging process, whereas the matrix is not invertible [3]. In the case of a vanishing standard deviation (), the application of the diagonal loading technique [3] guarantees the existence of the inverse data covariance matrix, since the condition number of is close to unity at the optimal diagonal loading parameter.
Furthermore, there exists a unitary transformation so thatwhereis a diagonal matrix that contains the eigenvalues of . Therefore, the inverse data covariance matrix results inwhereFurther transformation deliversso that the square root of the inverse data covariance matrix is defined aswhere
Statements
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.
Acknowledgments
The authors are grateful for stimulating discussions and helpful suggestions by Karl-Heinz Glassmeier, Patrick Kolhey and Alexander Schwenke.
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.
References
1.
CaponJ. High-resolution Frequency-Wavenumber Spectrum Analysis. Proc IEEE (1969) 57:1408–18. 10.1109/PROC.1969.7278
2.
ToepferSNaritaYHeynerDMotschmannU. The Capon Method for Mercury's Magnetic Field Analysis. Front Phys (2020) 8:249. 10.3389/fphy.2020.00249
3.
ToepferSNaritaYHeynerDKolheyPMotschmannU. Mathematical Foundation of Capon's Method for Planetary Magnetic Field Analysis. Geosci Instrum Method Data Syst (2020) 9:471–81. 10.5194/gi-9-471-2020
4.
MotschmannUWoodwardTIGlassmeierKHSouthwoodDJPinçonJL. Wavelength and Direction Filtering by Magnetic Measurements at Satellite Arrays: Generalized Minimum Variance Analysis. J Geophys Res (1996) 101:4961–5. 10.1029/95JA03471
5.
SuWGuHNiJLiuG. Performance Analysis of MVDR Algorithm in the Presence of Amplitude and Phase Errors. IEEE Trans Antennas Propagat (2001) 49(12):1875–7. 10.1109/8.982472
6.
NaritaY. A Note on Capon’s Minimum Variance Projection for Multi-Spacecraft Data Analysis. Front Phys (2019) 7:121. 10.3398/fphy.2019.00008
7.
SahraouiFBelmontGRezeauLCornilleau-WehrlinNPinçonJLBaloghA. Anisotropic Turbulent Spectra in the Terrestrial Magnetosheath as Seen by the Cluster Spacecraft. Phys Rev Lett (2006) 96:7. 10.1103/PhysRevLett.96.075002
8.
HuangSYZhouMSahraouiFDengXHPangYYuanZGet alWave 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. 10.1029/2010JA015335
9.
NaritaYComişelHMotschmannU. Spatial Structure of Ion-Scale Plasma Turbulence. Front Phys (2014) 2:13. 10.3389/fphy.2014.00013
10.
RobertsOWLiXJeskaL. A Statistical Study of the Solar Wind Turbulence at Ion Kinetic Scales Using the K-Filtering Technique and Cluster Data. ApJ (2015) 802:1. 10.1088/0004-637X/802/1/2
11.
WangTCaoJFuHMengXDunlopM. Compressible Turbulence with Slow-Mode Waves Observed in the Bursty Bulk Flow of Plasma Sheet. Geophys Res Lett (2016) 43(5):1854–61. 10.1002/2016GL068147
12.
ZhangLHeJNaritaYFengX. 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. 10.1029/2019JA027413
13.
NaritaYPlaschkeFMagnesWFischerDSchmidD. Error Estimate for Fluxgate Magnetometer In-Flight Calibration on a Spinning Spacecraft. Geosci Instrum Method Data Syst (2021) 10:13–24. 10.5194/gi-10-13-202
14.
WernerD. Funktionalanalysis. Berlin: Springer (2007)
15.
BelsleyDAKuhEWelschRE. 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, editorsGaussCFWeberW. Leipzig: Weidmannsche Buchhandlung (1839). p. 1–57.
17.
GlassmeierK-HTsurutaniBT. Carl Friedrich Gauss—General Theory of Terrestrial Magnetism—A Revised Translation of the German Text. Hist Geo Space Sci (2014) 5:11–62. 10.5194/hgss-5-11-2014
18.
ToepferSNaritaYGlassmeierK-HHeynerDKolheyPMotschmannUet alThe Mie Representation for Mercury's Magnetic Field. Earth Planets Space (2021) 73:65. 10.1186/s40623-021-01386-4
19.
DodelsonS. Modern Cosmology. London: Academic Press (2003).
Summary
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
Volume
9 - 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
Updates
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
This article was submitted to Space Physics, a section of the journal Frontiers in Physics
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.