Front. Phys.Frontiers in PhysicsFront. Phys.2296-424XFrontiers Media S.A.10.3389/fphy.2019.00008PhysicsPerspectiveA Note on Capon's Minimum Variance Projection for Multi-Spacecraft Data AnalysisNaritaYasuhito12*1Space Research Institute, Austrian Academy of Sciences, Graz, Austria2Institut für Geophysik und Extraterrestrische Physik, Technische Universität Braunschweig, Braunschweig, Germany
Edited by: Hermann Lühr, Helmholtz Center Potsdam German Geophysical Research Center (GFZ), Germany
Reviewed by: Peter Haesung Yoon, University of Maryland, College Park, United States; Xochitl Blanco-Cano, National Autonomous University of Mexico, Mexico
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.
Capon's minimum variance projection for the multi-point measurements is revisited using the method of likelihood function to derive the minimum variance projection and a simplified error estimate analytically. Theoretical construction of the minimum variance projection assumes a Gaussian form of the likelihood function and also regards the data covariance as a proxy of the noise covariance. The minimum variance projection is extended to the problem of two-spacecraft mode decomposition in the Mercury magnetosphere in which the magnetic field is a superposition of the constant field from the current sheet and the dipolar field from the planet. The extension of the Capon estimator (the data-variance projection) can identify the signal amplitudes of the different fields with a sufficient accuracy when the statistical averaging is properly done. The Capon estimator serves as a powerful analysis tool when the spatial resolution is limited to only a few points.
adaptive filter theoryCapon estimatormulti-spacecraft data analysiswaves and turbulencemode decomposition1. Capon's Minimum Variance Projection
Minimum variance projection introduced by Capon [1] has a wide range of applications in the geophysical and space physical research fields whenever multi-point measurements are available.
Various formations are possible in the multi-point measurements (Figure 1): the THEMIS mission in a one-dimensional array aligned with a magnetic field line [2], the Swarm mission spanning a plane with three spacecraft [3, 4], and the Cluster mission [5] the MMS mission [6] forming a tetrahedron.
Multi-point spacecraft configurations.
The projection works with various kinds of shape vectors (or models for the data) by minimizing the projection error without changing the amplitude of the signal amplitude. The minimum variance projection is, after Capon [1] or Haykin [7], obtained by imposing a constrained optimization:
minimize w→Rw→ subject to w→h→=1,
or, by formulating into a variational problem using the variation operator δ[⋯ ] and the Lagrangian multiplier λ,
δ[w→Rw→-λ(w→h→-1)]=0.
Here w→ is the weight vector operating on the measurement covariance matrix R=〈d→d→〉, d→ the measurement data in a vectorial form, 〈⋯ 〉 the operation of ensemble averaging, h→ the shape vector. The problem (Equation 2) can analytically be solved. The solution is, blue by treating the shape vector as a complex-number vector,
w→=R-1h→h→†R-1h→.
See, for example, Haykin [7] for the derivation. The projected power (squared signal amplitude) is
P=w→†Rw→=[h→†R-1h→]-1.
A larger set of algorithms has so far been developed in the frame of direction-of-arrival estimation in the adaptive filter theory. Many algorithms are recently reviewed by Khmou et al. [8], including beamforming method, Bartlett method, Capon method, linear prediction method, maximum entropy method, Pisarenko harmonic decomposition, minimum norm, MUSIC algorithm, propagator method, and partial covariance matrix method. In the multi-spacecraft wave-field analysis, three methods are most relevant: beamforming method, Capon method, and MSR method [9].
The beamforming method can easily be implemented to the data analysis, but on the other hand provides a lower resolution of the signal-to-noise ratio in the spectral analysis compared to the Capon method when only few data points are used in the analysis [10, 11]
The Capon method is versatile in the multi-spacecraft data analysis because the method can be extended in various ways, e.g., to the vectorial data set (the wave telescope technique) [12, 13] to different field types (e.g., electric field and magnetic field used in the k-filtering technique) [14], and to the mode decomposition (this paper).
The MSR method [11] uses both the Capon (or the wave telescope) method and the eigenvector-based method (the extended MUSIC algorithm), and provides an improved signal-to-noise ratio in the wavevector spectrum. The MSR method is optimized in the estimate of the total fluctuation energy (i.e., trace of the spectral density matrix) but not to the matrix elements. Also, the MSR method is constrained to the isotropic noise assumption.
So far, three shape vectors for the Capon method have successfully been applied to the multi-spacecraft data analysis blue in the field of space physics
Plane waves [1, 12, 14]:
hi(kx,ky,kz)=exp[i(kxrxi+kyryi+kzrzi)],
where {kx, ky, kz} are the three components of the wavevectors, {rxi, ryi, rzi} the spatial coordinates of the i-th sensor, and i the imaginary number unit.
Spherical waves [15, 16]:
hi(k,rxc,ryc,rzc)=[∑i=1n1|r→i-r→c|2]-1/2exp[ik|r→i-r→c|]|r→i-r→c|,
where k is the wavenumber, and r→c=[rxc,ryc,rzc] the coordinate of the spherical wave center. The summation in the normalization constant runs over the number of sensors n.
Phase-shifted waves [17]:
hi(Φ,rxc,Δx,ky)=exp[i(Φ arctan(rxi-rxcΔx)+kyryi)],
where Φ is the amount of phase jump, rxc the x-coordinate of the phase jump center, Δx the x-width (or range in x) of the phase jump around the center, and ky the y-component of the wavenumber.
Capon's minimum variance projection uses the measurement data to guide the projection by minimizing the estimated power during the projection (in spirit of minimizing uncertainty) yet keeping the gain. We address the question here, “Why is the inversion of the data covariance matrix R−1 used in the Capon projection (Equation 3) and in the spectrum (Equation 4)?” In this paper, we offer an answer to this question, that is, the Capon method uses the measured data as a proxy of the noise property upon the optimization procedure. The method of the likelihood function is introduced to give an alternative and more instinctive derivation of Capon's method. Moreover, by revisiting Capon's minimum variance projection through the likelihood function method, it becomes clear that one may project the measurement data not only onto a single shape vector but also onto a multitude of shape vectors, which will enhance the capability of the multi-spacecraft data analysis.
It is worthwhile to note that the maximum likelihood and the minimum variance aspects of the Capon estimator are already covered in detail in the original paper by Capon et al. [18] for a seismic array problem. In this paper, the essence of the Capon estimator is reviewed (in section 2) and the estimator is extended to the problem of mode decomposition. As an application, a new method is constructed (in section 3) for two-spacecraft measurements in the magnetosphere to identify the magnetic field of the current sheet origin and the dipolar field from the planet, which is relevant to the BepiColombo mission.
2. A View From Likelihood Function2.1. Scalar Field
Minimum variance projection can be formulated using the likelihood function L as follows (see also [19]). Consider a model for the observational data as
di=his+ηi,
where di is the data element at the i-th sensor (i = {1, 2, ⋯ , n}, so n is the number of sensors), hi the shape vector given a priori as a model, s the signal (in which we are interested), and ηi the noise at the i-th sensor. The noise is characterized by the n × n noise matrix N
Nij=〈ηiηj〉.
Again, the angular bracket 〈⋯ 〉 denotes the ensemble averaging over different realizations. The averaging is important because otherwise the determinant of the matrix vanishes and the matrix cannot be inverted.
We minimize the squared deviation between the model and the data, i.e., we minimize χ2 constructed (Figure 2) because of the steeper gradient toward the minimum (or the extremum) than that of the Gaussian curve as follows,
χ2=〈∑j=1n∑i=1n(di-his)Nij-1(dj-hjs)〉
with respect to the signal s. We now assume a Gaussian shape for the likelihood function,
L∝exp[-χ22].
The likelihood function L represents the probability of finding the true signal value for s under the given data set di.
Likelihood function on linear and logarithmic scales. Likelihood function is approximated to a Gaussian curve around the maximum (or the minimum of χ2=-lnL in the parameter space k.
Our goal is to estimate the signal amplitude (denoted as s~) for the shape vector for the given data set by finding a maximum of the likelihood function or equivalently by minimizing the deviation χ2.
We differentiate χ2 with respect to the signal s,
∂χ2∂s=-2〈∑i=1n∑j=1nhiNij-1(dj-hjs)〉.
Requiring that the derivative be zero (which corresponds to the extremum), ∂(χ2)/∂s, we obtain
〈∑i=1n∑j=1nhiNij-1(dj-hjs)〉=0.
Equation (13) can be arranged into the following form,
〈∑i=1n∑j=1nhiNij-1hjs〉=〈∑i=1n∑j=1nhiNij-1dj〉.
Note that the part in the angular bracket on the left hand side of Equation (14) is already statistically evaluated, so one may leave the angular bracket out in the following discussion. For a scalar quantity of s, one may obtain the estimator s~ as
In general, the data covariance matrix R=〈d→d→〉 and the noise covariance matrix N must be evaluated separately, or the noise covariance matrix needs to be known a priori to determine the signal covariance CS. Capon's minimum variance projection is obtained by imposing (or using) the data covariance matrix as the noise covariance matrix, N → R. An explicit calculation yields
CR=CNh→R-1h→CN=(h→R-1h→)-1h→R-1h→(h→R-1h→)-1=(h→R-1h→)-12.2. Vector Field
The method with the likelihood function can be extended to the vector field treatment in a straightforward fashion. We construct the data model as follows.
d→=Hs→+η.
The data are arranged into a long vector, d→=[dx1,dy1,dz1,dx2,dy2,dz2,⋯dzN]. The signal amplitude vector is s→=[sx,sy,sz]. H is the shape-and-pointing matrix.
H=[hx1000hy1000hz1hx2000hy2000hz2⋮00hzN].
η is the noise with the same construction as that of d→ (with 3 × n elements). The essential difference from the scalar-field treatment is that the covariance becomes a matrix, e.g., CS=〈s→s→〉 for the signal covariance matrix. Repeating the same procedure (by taking care of the vector and matrix operations), the signal estimator for the m-th component (k = {x, y, z}) is obtained as
s~k=∑ℓx,y,z∑i=1n∑j=1n(CN)ℓkHiℓNij-1dj,
where the summation on ℓ runs over the x, y, and z components, and that on i and j over the number of sensors. The signal estimator is given in a matrix notation as
s→~=CN〈HN-1d→〉
The noise covariance matrix CN is given by
(CN-1)ℓk=∑i=1n∑j=1nHiℓNij-1Hjk
or, in a matrix notation,
CN=[HN-1H]-1.
The signal covariance matrix, when using the data covariance matrix again as the noise covariance matrix, is
The minimum variance projection can be extended to multiple shape vectors. By doing so, it is possible to decompose the measurement data set into a spectrum of m different modes or shapes. We construct a model for the measurement data (scalar field) with a multitude of modes and a noise.
di=∑α=1mhiαsα+ηi,
where di is the measured field at the i-th sensor, m the number of the modes introduced into the model, h→α={h→(1),h→(2),⋯,h→(m)} is a set of shape vectors for modes α = {1, 2, ⋯ , m}, the symbol sα is the signal amplitude at each mode, and ni the noise at the i-th sensor.
Derivation of the minimum variance estimator is essentially the same as that for the single mode χ2 minimization. The estimator for the signal amplitude at the α-th mode is obtained as
s~α=〈∑βm∑i=1n∑j=1n(CN)αβhiβNij-1dj〉.
The noise covariance matrix is a projection of the inverse noise matrix N−1 onto different modes (hα and hβ),
(CN-1)αβ=∑i=1n∑j=1nhi(α)Nij-1hjβ.
The signal covariance matrix for the Capon-type projection is obtained as
(C)αβ=CNh→αN-1RN-1h→βCN→[hα→R-1h→β]-1.
The diagonal elements in the signal covariance matrix are the power (squared signal amplitude) of each mode.
2.4. Error Estimate
A useful form for the one-sigma error (68 % confidence) for the minimum variance estimator can be evaluated from the likelihood function. Ideally, the inverse of the noise matrix N−1 must be known, but the error estimate requires the knowledge on the noise property. Still, we obtain an insight by considering a special case, that is, we model the noise matrix N as diagonal with a value of each diagonal element CS+CN (which is a scalar). The likelihood function for a simplified error estimate is again modeled as Gaussian (see, e.g., Equation 11.21 in Dodelson [19]),
L∝1(CS+CN)n/2exp[-12〈∑i=1n(di-his)2〉CS+CN]
One-sigma error is obtained as the second-order derivative of the logarithm of the likelihood function (which is essentially χ2 in the Gaussian likelihood model):
σS=[-∂2∂CS2(lnL)]-1/2=2n(CS+CN)→8CS2n.
Here Equation (40) is evaluated at the peak of the likelihood function to obtain (Equation 41), and the signal covariance CS is used as a proxy of the noise covariance CN for Capon's minimum variance projection in deriving (Equation 42). Thus, the simplified error of the Capon-estimated signal power is σ2=8CS2/n. See Appendix for derivation of Equation (40).
3. Application
Mode decomposition using the minimum variance projection serves as a powerful analysis tool when the measurements are limited to only few spatial points. A test using a synthetic data set is presented with two-spacecraft measurements in the Mercury magnetosphere in view of the BepiColombo mission [20].
3.1. Setup
Magnetic fields are modeled as a superposition of two different fields (or modes), B(a) and B(b) and noise η in the magnetosphere at two sensor locations, r1 = 480 km (planetary orbiter) and r2 = 590 km (magnetospheric orbiter) above surface (at a radius of Rs = 2, 440 km for the planetary center). The measurement data are thus
where the first mode is the magnetic field from the current sheet
B(a)=s(a)=const.
and the second mode is the dipolar field from the planet
B(b)=s(b)(rRs)-3.
Measurements are assumed to be on the magnetic equatorial plane such that the magnetic fields have only one non-vanishing component (say, the z component). The geometrical configuration is illustrated in Figure 3. The true signal amplitudes are s(a) = 20 nT and s(b) = 200 nT, respectively. Noise is blue assumed to be Gaussian distributed with a standard deviation of σn = 1 nT. The goal of the numerical test using the synthetic data is to estimate the signal amplitudes for the two modes, s~(a) and s~(b) using the noise-variance minimization and the data-variance minimization.
Superposition of a constant field from the current sheet (mode 1) and a dipolar field from the planet (mode 2) in a region between the distance to the current sheet (from the planetary center) Rc and the that to the planetary surface Rs.
3.2. Analysis – Preparation
The shape vector for the first mode (a constant field) is
h→1=[11]
and that for the second mode (decaying field) is
h→2=[(r1Rs)-3(r2Rs)-3].
The shape matrix is constructed as [h→1h→2],
h=[1(r1Rs)-31(r2Rs)-3].
The data matrices are averaged over different realizations or samples,
R=〈d→d→〉,
using an averaging size of Ns. The measurement data vector is averaged when using the minimum variance estimator as 〈d→〉. We study the minimum variance estimators for different sampling sizes Ns. The measurement noise matrix follows uncorrelated Gaussian statistics,
N=〈η→η→〉=[σn200σn2].3.3. Analysis – Projection
The noise minimum variance estimator for each of the modes α = {a, b} is
(s→~(nv))α=(CN)αβ(h)βi(N-1)ij〈d→j〉.
The data minimum variance estimator for each of the modes α = {a, b} is obtained by replacing the measurement noise covariance N by the measurement data covariance R=〈d→d→〉 and also replacing the noise projection CN by the data projection CR as
(s→~(dv))α=(CR)αβ(h)βi(R-1)ij〈d→j〉.
The noise covariance matrix and the data covariance matrix are essentially the Capon projection, namely,
CN=[hN-1h]-1CR=[hR-1h]-1.
The noise minimum variance estimator needs the knowledge on the noise property and does not require the measurement data themselves. Therefore, the noise minimum variance estimator can be applied without the presence of the data and may be useful when planning a measurement of an experiment. For an uncorrelated Gaussian noise statistics, the noise minimum variance estimator computes the mode amplitude as a linear combination. The data minimum variance estimator, in contrast, requires the data but not the knowledge on the noise property. The mode amplitude is computed in a non-linear fashion, i.e., the measurement weight for each sensor is influenced by the data.
3.4. Results
Signal amplitudes for the two modes are obtained using the noise variance projection (Equation 52) and the data variance projection (Equation 52), and are graphically displayed as a function of the averaging size (or the number of realizations or samples) in Figures 4, 5 together with the one-sigma errors. The both estimators can find the true signal amplitudes (20 nT for the mode 1 and 200 nT for the mode 2) within the error bar, and the analysis using a larger statistical sampling size is beneficial in reducing the error. Yet, the accuracy is by far improved in the data-variance projection. The error bar is about 10 nT for an averaging size of 10–100 and becomes only a few nT or better (cf. the noise amplitude is 1 nT) for an even larger averaging size in the data-variance projection. The estimated amplitudes sufficiently converge to the true values for averaging sizes above 10, too. In contrast, the noise-variance projection still exhibits random deviations of the signal amplitudes from the true values for a larger averaging size (for example, size of 1,000).
Amplitudes for the mode 1 (constant field from the current sheet) and the mode 2 (dipolar field at the planetary surface) estimated by the noise-variance projection as a function of the averaging size.
Amplitudes for the mode 1 and the mode 2 (the same format as Figure 4) estimated by the data-variance projection.
4. Outlook
Capon's projection is a useful tool when the noise property is unknown, and has a higher flexibility for various applications and extensions compared to the beamforming or the MSR methods. The extension of the minimum variance projection onto a multitude of shape vectors in the data opens the door to a decomposition method for the multi-point data. An application is presented for a decomposition of the multi-point data into a constant magnetic field and a dipolar field in view of the Mercury magnetosphere. It is comforting that even two-point measurements are capable of identifying the signal amplitudes when a sufficient amount of data is obtained for the proper averaging operation.
Another application is a decomposition into a set of orthogonal function basis. Capon estimator and its extension to the mode decomposition can be applied to various wave fields in the solar wind, the foreshock, the magnetosheath, and the magnetotail regions as well as various static fields as far as the spatial structure can be properly modeled (such as the dayside magnetosphere as presented in this paper). For example, the source locator uses the lowest-order spherical Bessel function j0(x) = sin(x)/x and the lowest-order Neumann function n0(x) = −j−1(x) = −cos(x)/x. The expansion into a series of spherical Bessel functions or cylindrical functions (presumably with a cutoff) is a possible application to the multi-point measurements, e.g., identification or reconstruction of the spherical propagation or the vortical shape or motion.
Capon's minimum variance projection is not limited to the search for wave propagations or mode decomposition into different sources of the magnetic fields, but the method can be applied to solitary, spatially-localized structures. A useful example may be the KdV (Korteweg-de Vries) soliton (or ion-acoustic solitons in the case of plasmas) characterized by the following shape vector, hi(A,D,c,x0)=Asech2[D(xi-ct+x0)] where c is the phase speed of the propagation, A the amplitude, and D the width of the soliton structure around the peak. The amplitude and the width are determined by the phase speed. For ion-acoustic solitons, the amplitude is given as A = 3c and the width D=c/2 [21].
Author Contributions
The author confirms being the sole contributor of this work and has approved it for publication.
Conflict of Interest Statement
The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Discussions with Tohru Hada, Uwe Motschmann, and Daniel Heyner are greatly acknowledged to improve the quality of the manuscript.
ReferencesCaponJ. High resolution frequency-wavenumber spectrum analysis. Proc IEEE (1969) 57:1408–18. 10.1109/PROC.1969.7278AngelopoulosV. The THEMIS mission. Space Sci Rev. (2008) 141:5–34. 10.1007/s11214-008-9336-1Friis-ChristensenELührHHulotG. Swarm: a constellation to study the Earth's magnetic field. Earth Planets Space (2006) 58:351–8. 10.1186/BF03351933Friis-ChristensenELührHKnudsenDHaagmansR. Swarm–an earth observation mission Investigating geospace. Adv Space Res. (2008) 41:210–6. 10.1016/j.asr.2006.10.008EscoubetCPFehringerMGoldsteinM. Introduction to the cluster mission, Ann Geophys. (2001) 19:1197–200. 10.5194/angeo-19-1197-2001BurchJLMooreTETorbertRBGilesBL. Magnetospheric multiscale overview and science objectives. Space Sci Rev. (2016) 199:5–21. 10.1007/s11214-015-0164-9HaykinS. Adaptive Filter Theory, 2nd. ed., Prentice Hall Information and System Science Series. Upper Saddle River, NJ: Prentice-Hall Inc. (1991).KhmouYSafiSFrikelM. Comparative study between several direction of arrival estimation methods. J Telecomm Inform Tech. (2014) 2014:41–8.NaritaYGlassmeierKHMotschmannU. Wave vector analysis methods using multi-point measurements. Nonlin Processes Geophys. (2010) 17:383–94. 10.5194/npg-17-383-2010MotschmannUWoodwardTIGlassmeierKHDunlopMW. Array signal processing techniques. In: GlassmeierKHMotschmannUSchmidtR. editors. Proceedings of the Cluster Workshop on Data Analysis Tools, Braunschweig (1995). ESA SP-371, p. 79–86. p. 28–3.NaritaYGlassmeierKHMotschmannU. High-resolution wave number spectrum using multi-point measurements in space–the multi-point signal resonator (MSR) technique. Ann Geophys. (2011) 29:351–60. 10.5194/angeo-29-351-2011MotschmannUWoodwardTIGlassmeierKHSouthwoodDJPinçonJ. Wavelength and direction filtering by magnetic measurements at satellite arrays: generalized minimum variance analysis. J Geophys Res. (1996) 101:4961–6. 10.1029/95JA03471GlassmeierKHMotschmannUDunlopMBaloghAAcuñaMHCarrC. Cluster as a wave telescope–first results from the fluxgate magnetometer. Ann Geophys. (2001) 19:1439–47. 10.5194/angeo-19-1439-2001PinçonJLLefeuvreF. Local characterization of homogeneous turbulence in a space plasma from simultaneous measurements of field components at several points in space. J Geophys Res. (1991) 96:1789–802. 10.1029/90JA02183ConstantinescuODGlassmeierKHMotschmannUTreumannRAFornaçonKHFränzM. Plasma wave source location using CLUSTER as a spherical wave telescope. J Geophys Res Space Phys. (2006) 111:A09221. 10.1029/2005JA011550ConstantinescuODGlassmeierKHDécréauPMEFränzMFornaçonKH. Low frequency wave sources in the outer magnetosphere, magnetosheath, and near Earth solar wind. Ann Geophys. (2007) 25:2217–28. 10.5194/angeo-25-2217-2007PlaschkeFGlassmeierKHConstantinescuODMannIRMillingDKMotschmannU. Statistical analysis of ground based magnetic field measurements with the field line resonance detector. Ann Geophys. (2008) 26:3477–89. 10.5194/angeo-26-3477-2008CaponJGreenfieldRJKolkerRJ. Multidimensional maximul-likelihood processing of a large aperture seismic array. Proc IEEE (1967) 55:192–211. 10.1109/PROC.1967.5439DodelsonS. Modern Cosmology. London: Academic Press (2003).BenkhoffJvanCasteren JHayakawaHFujimotoMLaaksoHNovaraM. BepiColombo–comprehensive exploration of mercury: mission overview and Science goals. Planet Space Sci. (2010) 58:2–20. 10.1016/j.pss.2009.09.020IchikawaYHWatanabeS. Solitons, envelope solitons in collisionless plasmas. J Phys Colloques (1977) 38:15–26. 10.1051/jphyscol:1977603Appendix: Derivative Calculation
Calculation for Equation 40 is as follows. The first-order derivative of the Gaussian error likelihood function (Equation 39) with respect to the variance CS is obtained: