Post-processing scheme for modeling the lithospheric magnetic ﬁeld

. We investigated how the noise in satellite magnetic data affects magnetic lithospheric ﬁeld models derived from these data in the special case where this noise is correlated along satellite orbit tracks. For this we describe the satellite data noise as a perturbation magnetic ﬁeld scaled independently for each orbit, where the scaling factor is a random variable, normally distributed with zero mean. Under this assumption, we have been able to derive a model for errors in lithospheric 5 models generated by the correlated satellite data noise. Unless the perturbation ﬁeld is known, estimating the noise in the lithospheric ﬁeld model is a non-linear inverse problem. We therefore proposed an iterative post-processing technique to estimate both the lithospheric ﬁeld model and its associated noise model. The technique has been successfully applied to derive a lithospheric ﬁeld model from CHAMP satellite data up to spherical harmonic degree 120. The model is in agreement 10 with other existing models. The technique can be in principle extended for all kind of potential ﬁeld data with ”along track” correlated errors.


Introduction
All geophysical data are contaminated by signals that cannot be easily described by models.These poorly parameterized contributions are often treated as errors and they most of the time exceed the pure instrumental noise.These kind of errors are particularly difficult to deal with because they are often correlated in space and/or time.Further they may not follow a gaussian distribution.Yet, handling properly the data errors is at the heart of the data interpretation process and it usually requires their full statistical description -i.e. for a set of discrete measurements, the knowledge of the full covariance matrix of the data errors.
Geopotential data -i.e.gravity and magnetic measurements -are not an exception.For these types of data, the inverse problem that consists in finding the sources of the signals is particularly ill-posed, and the proper statistical description of the data errors is necessary.Failing to do so may lead to false conclusions about the signal sources.From a practical point of view, scientists have been relatively successful in estimating a priori the noise in gravity or magnetic data sets, however correlations between errors have been most of the time ignored.This is partly because, when known, the full covariance matrix for the data errors is generally so large that it cannot be handled easily, even on modern computers (but see Langel et al. (1989); Holme and Bloxham (1996); Rygaard-Hjalsted et al. (1997); Holme (2000) as examples where correlated errors are accounted for in geomagnetism).
The effects of these correlation errors are obvious in airborne, marine and satellite data.Typically, in all these type of surveys, the data are collected along linear paths and, after processing, the correlation errors become apparent as offsets between adjacent tracks.They then appear in maps and models as spurious anomalies, elongated in the direction of the tracks.An example of such an effect is shown in this manuscript for magnetic models derived from satellite data.The traditional way of dealing with this noise has been to perform a "leveling" of the data.In airborne geophysics, the approach mainly consists in deriving for each track a polynomial expression that is subtracted from the data such as to minimize data differences at the cross-over points (see for a review e.g.Hamoudi et al. (2010) ).The method has been also adapted to satellite magnetic data.In that case a large-scale field of external origin is fitted to a data set made of only few tracks.This allows to successfully derived magnetic field models of the lithosphere to relatively high degree.A well known example is the MF series of models -e.g.Maus et al. (2008).However the method, as applied to satellite data, has its drawbacks.The effects of its application have been carefully studied in Thebault et al. (2012) and it appears that, depending on the way the method is applied, it can lead to significant distortions of the final model.However, the weakest point of this so called "along track filtering" approach is the impossibility to estimate how much the processing applied distorts the model.For this aspect, post-processing techniques are preferable.
So far post-processing techniques have been developed and applied only to models derived from satellite gravity data -e.g.Kusche (2007).To the authors knowledge such techniques have never been applied to magnetic models although we should note the attempt to estimate the model covariance matrix in Lowes and Olsen (2004).In this manuscript we present and apply such a postprocessing scheme for a model of the magnetic lithospheric field derived from ten years of CHAMP satellite data (Reigber et al., 2005).Although we are presenting this work from its application side, it has deeper roots: We investigated how typical noise correlated patterns leak, through a least squares fitting process, inside a magnetic model of the lithospheric field.Therefore this leads to a model of the noise inside the lithospheric model.Once such a noise model is available, numerous postprocessing schemes are possible; we just applied one specific approach to show that the noise model we obtain is relevant.The final resulting model of the lithospheric field is nonetheless of high qual-ity and compares well with other recently released models (e.g.MF7 that is not published -but the MF6 is presented in Maus et al. (2008), CHAOS-4 (Olsen et al., 2010b)) as well as older models (see for a review Thébault et al. ( 2010)).
The manuscript is organized as follows.In the next section we set the hypothesis and approximations, derive the general expression for the noise model and give examples of possible noise depending on the characteristic of the perturbation magnetic field in the data.In the third section we describe in details the two steps process towards the final lithospheric field model, the resulting model is then discussed.We conclude in the last section.

The lithospheric noise model
In this section we present a noise model for a lithospheric model estimated from a set of radial magnetic data.We choose to present this case only in the main part of this manuscript as the equations are relatively simple to derive.The description for the usual case where the lithospheric model is obtained from the three components of a magnetic data set is given in appendix A.

Theory
We consider a magnetic data set made of radial component readings along a single CHAMP satellite half-orbit during night times.For simplicity we will assume that a track follows a meridian -i.e. it corresponds to a single longitude value , which is a reasonable approximation for near-polar orbiting satellites.Several magnetic field sources are contributing to these data (Hulot et al., 2007), typically the core and lithospheric fields, the ionospheric and magnetospheric fields, and the fields generated by field aligned currents.Other contributions exist, as the field induced in the conductive layer of the Earth, but they are of much weaker amplitudes.Mathematical models are available for all these contributions and can be subtracted from the data, leaving residuals due mainly to the limited precision of these models.In particular, the description of the external field is not very accurate and the residuals obtained along that half-orbit track contain relatively long wavelengths.We assume that these residuals are well approximated by the radial component of an external magnetic field model that does not present time dependencies.It is hereafter named as the perturbation field and writes: where ϵ k n is the Gauss coefficient of degree n and order k, a = 6371.2km is the Earth's reference radius, Y k n (θ,ϕ) are the Schmidt semi-normalized Spherical Harmonics (SHs).We use along this manuscript the convention that negative orders, k < 0, are associated with sin(|k|ϕ) terms whereas null or positive orders, k ≥ 0, are associated with cos(kϕ) terms.
We consider also a model of the radial component of a magnetic field of internal origin with no temporal dependencies.This model becomes below the lithospheric noise model we want to derive: (2) It is not possible to separate external field contributions from internal field contributions for data collected along a single meridian (Olsen et al., 2010a) -i.e. a single half-orbit, hence we can fit by least-squares the residuals defined in Eq. 1 with the lithospheric model given in Eq. 2 and find a non-zero solution.This least-squares solution is found by minimizing the functional: where θ i are sampling points along the half-orbit, ϕ j is the longitude of the meridian that we labeled with the subscript j and w i are weights that are defined below.
Over 10 years, the CHAMP satellite has collected data along a large number M of half-orbits.We assume now that for each orbit the perturbation field model defined by Eq. 1 is scaled by a number η j and that all orbits are at the same radius r.This latter point is clearly a strong approximation but there is no obvious way to avoid it.Again, these external field contributions can be interpreted as a field of internal origin.To estimate this field, the functional we have to minimize is then: Minimizing Φ for the Gauss coefficients gm l leads to a system of equations: where g = [g m l ] {l,m} .The matrix product A t A is derived from Eqs. 2, 4 and the elements of this product associated with the degrees and orders l,l ′ ,m,m ′ writes: where the product The variable Π mm ′ has been introduced to cover three cases: and is symmetric relative to its subscripts -i.e.
The elements of the right hand side vector of Eq. 5 are: Depending on the sign of the orders m and k, χ k m takes the following values: and, as for Π mm ′ , its is symmetric relative to its subscripts: χ k m = χ m k .We note at this point that it is important to have the η i constant along half-orbits otherwise the summations over latitudes and longitudes could not be separated.
For very large number M of orbits uniformly distributed along longitudes, the quantity Π mm ′ tends to a δ-function -i.e.Π mm ′ ≃ ( 1 2 + 1 2 δ m0 )δ mm ′ .Further, by setting the weights w i to w i = sinθ i and assuming that the sampling points are evenly spaced over the full meridian, we have 2l+1 δ ll ′ .The value given to the weights w i is less important than insuring the orthogonality of the legendre functions through the product ⟨P |m| l ,P |m| l ′ ⟩.This is also what a modeler tries to acheive when building a lithospheric magnetic field model from real data.However, assuming that both Π mm ′ and ⟨P |m| l ,P |m| l ′ ⟩ can be regarded as δ-functions, it is easy to see from Eq. 6 that the product matrix A t A is diagonal.Now turning to Eqs. 8, 9, if the η i form a set of uncorrelated random variables, the χ k m are also random variables with zero mean.The Gauss coefficients for the lithospheric noise model in Eq. 2 are then obtained by combining Eqs. 5, 6 and 8.They are: They correspond to the noise in a lithospheric field model that would be generated by un-modeled external fields in the radial component of magnetic data.Similarly, it is straightforward to find the noise in a lithospheric field model (i.e.static internal field model) generated by a perturbation field of internal origin.This case is relevant for signals generated in the lower E-region ionosphere (e.g. at 110 km altitude) when data are acquired at satellite altitudes.Other possible sources for this type of noise are the un-modeled induced fields generated in the conductive layers of the Earth by rapid variations of the external fields.It gives: where ı k n are the Gauss coefficients for the ionospheric and/or induced field models.In order to understand the behaviour of the lithospheric noise model, it is important to have an estimate of the probability density function of the random variable χ k m .Assuming the random variable η is normally distributed with variance v η then χ k m appears to be also normally distributed.The set of χ k m are uncorrelated with the exception that χ k m = χ m k .Further the χ k m have a variance v χ that depends on v η , the number of half-orbits M , the orders k and m.Possible values of the variance v χ are given in table 1.These variances have been derived from numerical experiments involving 20000 independent realizations of the random variables χ k m calculated from the same number of uniformly distributed orbits.They can also be estimated analytically as shown in appendix B. Figure 1 presents the histograms for few values of m and k.
In the remaining parts of this section we consider only the noise model given by Eq. 11.The general behaviour of the noise characterized by Eqs.11 and 10 is basically the same.In particular they have the same dependence relative to the degree l.These two noise models are only relevant for the cases where the radial components of vector data are used.The way the noise propagates in a lithospheric model is different if the three components of the vector data are fitted.The corresponding equations for that case are relatively complex and given in appendix A.

Examples
In order to understand the main characteristics of the noise model defined by Eqs. 2 and 11, we present in this section the results of forward modelling calculations for a given choice of Gauss ⟩ are calculated numerically.These products are relatively difficult to estimate accurately as the P m l (x) functions are oscillatory.However, an adaptive Gaussian quadrature (Piessens et al., 1983;Kahaner et al., 1989) was ultimately chosen as it gave the best results.

Dipole perturbation field
For this first example we use a simple model for the perturbation field of internal origin made of a single spherical harmonic n = 1, k = 1.Specifically, we set ı 1 1 = 1 and ı k n = 0 for {n,k} ̸ = {1,1}.This type of noise in satellite data could result from a poor modelling of the field induced by a large-scale external field in the conductive layers of the Earth.In that case Eq. 11 reduces to: and the noise in the radial component of the field of internal origin is: where r ′ is the modelling radius that is set to r ′ = a = 6371.2km in this example.As the observation radius r is expected to be larger than the modelling radius, the short wavelengths dominate the model due to the ratio r r ′ raised to the power l − 1 in the right hand side of Eq. 13.
In Fig. 2, the model defined by Eq. 13 is mapped for the model coefficient ı 1 1 = 1 nT, an observation radius at 300 km altitude (r = 6671.2km) and the random variables χ 1 m with variances defined in table 1 using v η = M .The maximum SH degree involved is L = 120.We observe that the noise model is symmetric relative to the equator, vanishes at the poles, and is made of East-West oscillating anomalies typical of the noise in lithospheric field model drived from satellite data.We note that these characteristics are independent from the sign of the SH order k as only the random variable χ k m depends of this sign in Eq. 13.The obtained symmetry of the model is due to the product ⟨P The power spectrum of the model calculated at r ′ = 6371.2km is also plotted in Fig. 2. It presents some variability due to the use of a single SH in Eq. 12. Nonetheless, the behaviour is generally along a ( r r ′ ) 2l trend as it would be expected for a white noise at satellite altitude.Although the small wavelengths overshadow the larger wavelengths, the latter are also present in the noise model.It is clear that any magnetic field model derived from satellite data is contaminated by such a noise at all wavelengths unless pertinent processing steps are applied.

Auroral electrojet and field aligned currents
Another expected source of noise in satellite data is associated with auroral electrojet and/or associated field-aligned currents.We do not aim at a precise description of the disturbance field but just consider the radial component of a perturbation field of internal origin, mapped in Fig. 3 (left), and defined by : The model was built simply by defining a circular ridge in a geomagnetic system of coordinates, that was then rotated in the usual reference frame.We recall that in our approach this field is scaled by a random variable with zero mean for each orbit.Therefore it is more the geometry of the field that is important here than its true value.We see that the perturbation field model is centred on the geomagnetic North pole and takes relatively large values up to 60 o colatitudes.The lithospheric noise model we obtain, derived from Eq. 11, is mapped in Fig. 3 (right).This model is also fairly well localized in latitudes as it basically vanishes in the southern hemisphere.However, it seems that the noise is propagating over all longitudes.The power spectrum of the model has essentially the same characteristic than in the previous example.The results of this example have to be analysed with some caution since real satellite orbits deviate from the exact polar direction at high latitudes.Nonetheless, we take out from these results that there is no need to describe precisely the longitudinal dependence of the rapidly varying field to obtain a realistic noise model.Therefore, in Eq. 11, the range of SH order k can be restricted to small values -e.g.k max = 2, even if the maximum SH degree in the model remains large -e.g.N = 30.This will reduce even further the number of parameters needed to describe the noise model.
3 Application to magnetic models of the lithosphere We are now using the results presented in the previous section to derive a model of the magnetic field generated in the lithosphere from real CHAMP satellite data.The process we applied to calculate such a model is in two stages.First we estimate a rough lithospheric field model from satellite data using an usual approach and a straightforward least-squares process (e.g.Lesur et al. (2008)).
Second, in the post-processing stage, we co-estimate a new lithospheric field model and a model of the noise where the output model of the first stage is used as data.The final results depend on the processes applied during the two stages and therefore both are described in independent subsections below.
In order to avoid confusion between the different models we use the following notations for the fields and Gauss coefficients: -the noisy lithospheric model, output of the first stage, is denoted using a .-e.g Bi for the magnetic field vector, -the lithospheric field model output of the post-processing (second) stage does not have any distinctive sign -e.g B i .
-the noise model is denoted, in the same way as in the previous section, using a .-e.g Bi .
3.1 First stage: Data set, data selection, model parameterization and model estimation Three component vector magnetic readings acquired during the ten years of the German CHAMP satellite mission are used.The data are selected for night-times and magnetically quiet days, in the same way as data are selected for the GRIMM series of core field model (Lesur et al., 2008(Lesur et al., , 2010)).
However, here the three components of the vector data are used and data in single star camera mode are rejected, whereas in the GRIMM selection scheme only the X and Y SM components are selected at mid and low latitudes.A core field model and a model of the large-scale external field with its internally induced counterpart are subtracted from these data, leaving mainly the contributions from the lithosphere and the noise.The core field model and external field models used are resulting from the derivation of GRIMM-3 (Lesur et al., 2011), but this is not seen as an important point in the processing: another core field model would have been possible -e.g.CHAOS-4 (Olsen et al., 2010b).
Next a first lithospheric field model up to SH degree 60 is derived, but our aim here is to reject outliers.The data corresponding to residuals larger than 3 times the standard deviation are rejected.
The value of the threshold, for each data type, is given in table 2. This selection process is known to potentially affect strongly the final lithospheric field model.At mid and low latitudes only few data are rejected, and those rejected data do not present clusters: No major difficulties are therefore expected there.At high latitudes however, a large amount of data are rejected and it is not possible to assess at this point if magnetic anomalies are erased or minimized there.We checked however, that outside the polar gaps due to the satellite orbits, the final data density at satellite altitude is everywhere large enough to allow for a lithospheric field model to be estimated up to SH degree 120.
In order to avoid spurious oscillations of the lithospheric model, further vertical down component data values were added over the polar gaps at an altitude of 6371.2 km.The data values were arbitrarily set to zero, and the associated weights for the inversion process were adjusted such that the lithospheric field model remained smooth while the misfit to the original satellite data stayed unchanged.We have tested other possible approaches, but the one we used gave the best results.
One could alternatively use vertical down component values derived from aeromagnetic maps.
The data set resulting from this selection process still consists in some 5 014 325 data values.A model of the lithosphere magnetic field, defined by Eq. 15 below, was fitted through a simple leastsquares process to the data.The data weights are set to mimic a homogeneous data repartition and therefore depend only on the inverse of the data density.
The power spectrum at the Earth surface of the resulting lithospheric field model is presented in The vertical down component of the model -i.e. the Z component, is mapped at 300 km above the Earth surface in Fig. 4, right.At this altitude the long wavelength lithospheric signal dominates but the noise is clearly visible, mainly over oceans, as elongated anomalies in the north-south direction -e.g. to the south of Australia.We point out that there are strong correlations between the estimated Gauss coefficients of the model and therefore the model cannot be truncated at an arbitrary degree without introducing artefacts.

Second stage: Model post-processing
The post-processing part consists in fitting a model of the magnetic field generated in the lithosphere together with the model of noise to a 300 km altitude map of the vertical down component of the field model Bi (θ,ϕ,r) (see Fig. 4).The noise model Bi we used is derived in appendix A and is parameterized by the variable χ k m and the Gauss coefficients of the perturbation model i k n .This inverse problem that consists in fitting the noise model and the lithospheric field model to BZi values presents some difficulties that are described first, results are given in a second subsection.

Inverse problem
We map the vertical down component of the magnetic field model Bi (θ,ϕ,r) at 29161 positions on a Gauss-Legendre grid at r = 300 km altitude.These data values are related to the Gauss coefficients g m l of the field model B i by the relation: where BZi (θ i ,ϕ i ,r) is the vertical down component of the noise model derived in appendix A, and ϵ i is an unknown noise.As the maximum SH degrees in Bi and B i are the same, it is clear that the g m l can be estimated such as B i fits exactly the values of BZi (θ i ,ϕ i ,r) with the noise model and the ϵ i not contributing to the problem.These latter contributions become necessary only when a priori smoothness requirements are introduced on B i .Hence the inverse problem consists in minimizing the functional Φ defined by: The first term insures the fit to the data BZi (θ i ,ϕ i ,r) whereas the second minimizes the integral of the squared horizontal gradient of the radial component of B i over a sphere of radius a =6371.2km.
The parameter λ controls the smoothness constraint applied on B i .
As stated above, the noise model Bi (Eqs.A3 and A12) is parameterized by the variable χ k m and the Gauss coefficients of the perturbation model i k n .A possibility is to set the perturbation model coefficients i k n , such that the model corresponds to a dipole field, and to try to estimate the χ k m .The inverse problem is then linear.However, for such a choice the derived lithospheric field model B i appears to be still contaminated by noise, probably because the perturbation models is more complex than a simple dipole.Therefore, there is no other option than to co-estimate the χ k m and i k n values.
As these quantities enter as products in Eq.A12, the inverse problem is non-linear and must be solved iteratively.We want to point out that finding the χ k m and i k n values in Eq.A12 or in Eq. 11 are two different problems with their own specific null-space and difficulties.In particular, i k n and i −k n values cannot be estimated independently if Eq. 11 is used.With Eq.A12 this estimation becomes possible solely because of the way the Y component data affect the noise model.However in both cases the maximum value for n can be relatively large, whereas the maximum value of k has to be small.We used in this work a maximum value of n: N = 20 and a maximum value for k : K = 1.As noted in section 2.2.2 most of the complexity in longitude of the noise model is carried by the χ k m ; there is no need for a large longitudinal complexity of the perturbation model.With such settings, the number of unknowns describing the noise model in Eq.A12 is reduced to N (2K +1) + K − K 2 for the i k n (i.e. 60 values) and (2K + 1)(2L + 1) − 2K 2 for the χ k m (i.e.721 values for L = 120).These numbers have to be compared with the number of unknowns in the lithospheric field model The iterative inversion process we followed to reach the solution presented in the next subsection is described in three steps: step-1 Find the g m l in Eq. 16 minimizing Φ (Eq.17) imposing χ k m = 0 for all possible m and k values.
step-2 Keeping the g m l unchanged, and starting with i k n = 1 for all possible n and k values, find iteratively the i k n and χ k m that minimize Φ in Eq. 17.
step-3 Iteratively find the g m l , i k n and χ k m that minimize Φ in Eq. 17, starting from the output of the step 2.

Results
The results were obtained by iteratively minimizing the functional defined in Eq. 17 following the process described above, with the parameter λ set to λ = 4.010 −5 such that the resulting field model has a power spectrum in its expected range.The level of noise is larger at high latitudes in the BZi (θ i ,ϕ i ,r), we therefore weight the data by 1 6 for magnetic latitudes higher than 50 ø .
The output of the step-1 described above, is a smoothed model obtained without co-estimation of the noise model.The map of this model vertical down component at radius 6371.2 km is shown in Fig. 5.The perturbation due to the along track noise in the satellite data are strong, particularly over Antarctica, and in the Indian, Atlantic and eastern Pacific oceans.This map is given here as reference for comparison with our final model obtained by co-estimation with the noise model.
The residuals to the fit to the data after the last step of the fitting process are mapped in Fig. 6, left.
The largest anomalies, as the Bangui anomaly in central Africa or the Kursk anomaly in western Russia, are clearly identifiable on this residual map, although they are not associated with too large residuals.There are very large clusters of residuals at high latitudes, and some of these residuals obviously correspond to lithospheric magnetic anomalies -e.g.North America, Southern tip of Greenland, Northern Europe.This is an incentive to work with localized system of representation and to define local constraints.Here we want to keep the processing as simple as possible and did not follow such approaches.It should be noted, however, that the amplitude of the residuals are clearly smaller than 1.5 nT and that there is only few traces of the "along track" noise in these residuals.
The effect of the smoothing on the model remains acceptable.
Figure 6, right, shows the power spectra of the field model B i and of the noise model Bi .Also plotted is the spectrum from MF7.The damping parameter λ in Eq. 17 has been adjusted to λ = 4.010 −5 such that the power spectrum does not present excessively high values at high degrees.
Overall, the derived map has the same level of energy than MF7 up to degree 100.Above that degree the spectra is clearly decreasing.Our opinion is that we are reaching at these SH degrees the maximum "global" resolution of the CHAMP data selected and processed following the technique described above.Improvements are probably still possible locally, particularly above the largest anomalies seen as Bangui and Kursk anomalies.The anomaly patterns are not as clearly defined as in MF7, but in numerous areas -e.g. the northern pacific, the map resulting from our processing is remarkably detailed.However, in the present case, the only difference with regards to a straightforward least-squares approach is the co-estimation of the noise model.In particular, there are no pre-processing steps such as data levelling (or microlevelling) with mostly unknown consequences on the final map, and the only data used are the CHAMP satellite data.We have made numerous experiments, and it appears that the determinant step for the final quality of the map is the data selection used to build the model Bi .Out of all these trial, the maps presenting the lowest level of noise are systematically the outputs of the step-2 of our processing.We decided not to show these results here because they are not consistent with the noise model presented in appendix A that assumes a model derived through a non-regularized scheme.It is however an approach worth studying.There are no major difficulties in estimating what the noise model should be for a lithospheric model built using a regularized least-squares process.

Conclusions
We have calculated the gauss coefficients describing the noise leaking in lithospheric magnetic field models when derived from satellite data.The noise models were derived to cover two cases: first when exclusively the radial components of the satellite data are used and second when all three components are used.The first case would be primarily applicable to gravity data, whereas the second, as we used it here, is better suited for magnetic data, although applications to vector gravimetry or gradiometry may be possible.We made several strong hypotheses to obtain these results.Particularly, we consider that the orbits are exactly polar, that they are at constant radius and that the sampling rate along an orbit is "ideal" -i.e. the relation We also make the assumption that the lithospheric field model is derived through a simple un-regularized least-squares process.This latter approximation is well verified for our application, but the three former are rather rough -e.g. in our data set the altitudes varies between 480 km and 250 km.However, it appears that final result does not suffer too much from these hypotheses.We insist here on the fact that the noise models do not represent the expected noise in the satellite data but the noise leaking in the derived lithospheric models.
It is interesting to notice that the amplitude of the noise generated depends on the variance of the random variable χ k m , that itself depends on the variance of the external field scaling factor η and the number of orbits M (see table 1).Therefore the usual choice of rejecting a significant part of the data because of its level of noise is questionable.For example when dealing with magnetic data, rejecting a full year of satellite data because of the high level of magnetic activity is unlikely to reduce the noise level in the model since the ratio v η M generally does not get smaller.We cannot comment however on a data rejection criteria based on the satellite altitude.
Another remarkable property of the noise models is their weak dependence with regard to the source of the noise.We used here perturbation models either from internal or external origin, but both lead to similar noise models.The same developments could be done for a noise described by spherical harmonics without reference to any specific source.For the case where only radial component data are used (Eqs. 10,11), such an hypothesis would not make any difference.
In the application to real data, the noise models were estimated in a post-processing scheme.The reason for this choice is that we did not know what kind of perturbation model B p (θ,ϕ,r) should be used.We have seen that for deriving a lithospheric field model, a dipole perturbation field is not leading to the best results.In an ideal case where the perturbation model is known, the best approach to the problem would be to build a covariance matrix C n for the noise from the variances given in table 1 and the Eq.A12.Such a covariance matrix can then be used as a regularization matrix in the least-squares fit of the lithospheric field model to the satellite data.However, even if the information provided by the estimated variances has not been used in our post-processing scheme, the resulting lithospheric field model is nonetheless much improved compared to what can be obtained through a simple smoothing (see the differences between Fig. 8 and Fig. 5).
One can question if parts of the lithospheric field can be removed by our post-processing steps and contribute to the noise model.The lithospheric noise model derived is only a combination of spherical harmonics with some strong correlations between the Gauss coefficients.Therefore, there is no doubt that part of the true lithospheric magnetic field model can contribute to the noise model.
It is however not possible to estimate a priori what this part is because it clearly depends on both the noisy lithospheric field model on which the post-processing is applied and the true lithospheric field we want to estimate.In order to test our scheme we have first applied the processing on a synthetic data set built on a Gauss-Legendre grid where both the lithospheric field model and the noise are known.We used only the radial component of the field and verified that the noisy lithospheric field model derived from these data was contaminated by a noise corresponding exactly to Eq. 11.However, the full inversion process revealed that part of the lithospheric field was seen as noise.We also applied the step-2 of our processing using a noise free synthetic lithospheric field model and the noise model defined by equation Eq.A12.Here again, despite a noise free data set, the lithospheric field is partly interpreted as noise.The amplitude of the obtained noise model is relatively large where there are strong, more or less aligned along orbits, anomalies.This is the case for the Kursk anomaly (51 • N, 37 • E), whereas the Bangui anomaly (4 • N, 16 • E) is apparently not affected by the processing.Outside few localized areas, the noise model remained relatively small.This impossibility to properly separate the noise from the lithospheric field is a common limitation of all existing processing methods.In our specific scheme, the only way we can reduce this effect is by constraining the perturbation model.We therefore recommend that the post-processing is used only when the noise is clearly identifiable at the smallest wavelengths, constraining this way the perturbation model.Overall, the proposed post-processing performs probably better than other smoothing techniques but that has to be tested on a case-by-case basis.
The work presented here opens numerous possibilities for processing data acquired along linear paths, as satellite data.The major difficulty when dealing with large data set, is to handle the correlated errors.Facing this problem we have here simply calculated how this correlated noise affects the derived model through a least-squares process.Extending this to regularized least-squares approaches is certainly possible.The same technique can be applied for calculating small-scale secular variations from satellite data, or to process yearly estimates of the core field.The technique is also applicable for airborne data using any local system of representation rather than spherical harmonics.Interesting developments are possible through the design of local filters.The link with oriented wavelets on the sphere is also promising.
With this definition it is clear that where v χ is the variance of χ k m and C N = v η I d , with I d being the identity matrix.We therefore obtain: and using the identity cosacosb = 1 2 (cos(a − b) + cos(a + b)) it follows: Over the CHAMP mission there is a large number of orbits and the ϕ i are uniformly distributed between [0;π].It follows that the last term on the right-end-side vanishes unless k = 0 or m = k = 0, and it can be verified that: as long as n is not too large.Therefore it comes:
coefficients ı k n and one realization of the set of random variables χ k m .The products ⟨P |k| n ,P |m| l

⟩
that vanishes if the Legendre function P |m| l is anti-symmetric -i.e.l − |m| is odd.An anti-symmetric model, vanishing at the equator but not at the poles, would have been obtained if ı 0 1 = 1nT would have been chosen in place of ı 1 1 = 1 nT.These symmetry/anti-symmetry characteristics are specific to models derived from the radial component alone.It can be seen in appendix A that these characteristics are lost when a noise model is obtained from the three vector components.

Fig. 4 ,
Fig.4, left, together with the power spectrum of the CHAOS-4 model.Both models present very similar spectra up to SH degree 60 or 65.Our model presents slightly less power around degree 70, possibly due to the selection technique used.Above SH degree 85, CHAOS-4 spectrum is strongly minimized, whereas our model presents a spectrum rising to high values, evidence of the predominance of noise in the model at these SH degree.The final misfits to the data are given in table 2.

Figure 7
Figure 7 maps, on the left, the noise model Bi , and, on the right, the perturbation model defined in Eq.A1.The noise model presents the expected East-West high frequency oscillations.The map cannot be directly compared with Fig. 5 because the patterns of the oscillations in Fig. 7 correspond to the noise present in Bi : the Fig. 5 is only a smoothed version of it.The perturbation model (Fig. 7, right) is dominated by a dipole term consistent with un-modelled contributions generated in 1-D conductive layers of the Earth by a large-scale, rapidly varying external field.Although this large-scale field is dominant, higher spherical harmonic contributions exist in the perturbation model and are determinant for the success of the post-processing.Our final result is a map of the vertical down component of the lithospheric field calculated at the Earth's surface (see Fig.8).The map includes all SH degrees of the lithospheric field model.The model is displayed with two different central meridians for a better view of the anomaly patterns.

Fig. 1 .
Fig. 1.Histograms of the random variable χ k m for several values of k and m.Is also plotted the dashed curve M •S √ 2πv χ exp{−e 2 /(2v χ )} where S is the histogram step length.

Fig. 4 .
Fig. 4. Left: Power spectra of the lithospheric field model (solid line) and of CHAOS-4b model (dashed line) calculated at the Earth's surface (i.e.r = 6371.2km).Right: Mapping of the vertical down component of the (noisy) lithospheric field Bi at r = 6671.2km.The largest magnetic anomalies dominate, but the "along track" noise is nonetheless visible over oceanic areas.

Fig. 5 .Fig. 6 .
Fig. 5. Map of the vertical down component of the lithosphere magnetic field model at r = 6371.2km radius derived after the step-1 of the processing chain.This corresponds to a smoothed model without co-estimation of the noise model.It is given here as a reference to be compared with figure 8. Along track noise is particularly visible around Antarctica, and in the Indian, Atlantic and eastern Pacific oceans.

Fig. 7 .
Fig. 7. Left: Map of the vertical down components of, left, the noise model, right, the perturbation model.Both maps have been calculated at r = 6371.2km radius.By definition the perturbation model is very smooth in longitude, but that does not preclude a large complexity for the noise model.

Fig. 8 .
Fig. 8. Map of the vertical down components of the final lithospheric field model B i .The map has been calculated at the Earth's surface (6371.2km).Although some noise is still visible in the northern Atlantic and over the southern polar cap, the noise level over mid latitudes has been greatly reduced.Anomalies are particularly well defined over continents, and Indian and Pacific oceans.

Table 1 .
Estimated variance of χ k m

Table 2 .
Thresholds and misfit values obtained when estimating Bi from selected satellite data.Mid and low latitudes are defined by magnetic latitudes in between ±55deg.Values are given in nT.