Abstract
This contribution deals with challenges encountered in real-world geomagnetic measurements and is focused on improving the performance of two variometer stations of Kelčany and Polom, which have been recently established in the Czech Republic. It is shown that a carefully designed full-field instrument, despite lacking temperature stabilization, can provide vectorial and scalar data accurate to a few nT, if the raw data were postprocessed by compensating for gain temperature coefficients—we show how this can be obtained by a precise calibration and long-term scalar measurements. We also show a method for suppressing nT-level spikes in the data due to nearby car traffic, by utilizing gradiometric measurement for detecting the car occurrences and by employing a linear optimization problem in order to find the parameters of the moving magnetic dipole and compensate for it. In this manner, we were able to reduce the anthropogenous noise due to car traffic while keeping as much original information as possible.
1. Introduction
Vectorial magnetometers which serve for monitoring of Earth’s field variations due to diurnal field changes, geomagnetic storms, etc., are standard instruments deployed at geomagnetic observatories and variation stations; they are mostly based on fluxgate sensors [1]. To achieve the best magnetometer (variometer) performance, it will usually be installed in a temperature-stabilized environment, either by the use of an active and magnetically clean heating-system or by selecting a highly temperature-stable location, preferably underground. This approach is of course demanding on the site selection and/or the necessary infrastructure. However, the variometer performance improved recently not only in terms of noise (the current state-of-the-art noise limit of fluxgate sensors is about 3–5 pT/√Hz @ 1 Hz) but also in temperature stability [1–3]. Also, with the advent of dc-precise 24-bit A/D converters, it is possible to build a “full-field” instrument which not only can monitor the magnetic field variations but also can provide the vector magnitude (scalar) from the three vector components. To achieve this, precise magnetometer calibrations are needed [4, 5], and the calibration parameters need to be long-term stable. The calculated total field value can be further used for temperature compensations.
An important aspect of real-world deployment of variometers is the anthropogenous noise at the selected site, which is at least in the Central European region difficult to obey by placing the instrument in a remote locality—due to extensive urban development, DC-railways, pipelines, etc. [6]. Thus, the anthropogenic noise should be estimated and, if a better location is not viable, a compensating or at least a detection method should be developed; the latter is the case mainly if the occurring disturbances are on local-scale, i.e., car traffic—we will show this is the case of one of our localities.
The following results were obtained from magnetometers running at three different localities in the Czech Republic. The reference low-noise data was obtained from the established INTERMAGNET geomagnetic observatory Budkov (BDV) in Southern Bohemia, which employs passive and active temperature stabilization. Two variometer stations were recently established—Polom (PLM) in Eastern Bohemia at the Czech-Polish border and Kelčany (KEL) in south Moravia—see Figure 1. The motivation to have all three stations is clear—first, having data redundancy is important for magnetic field observations, forecasts, and data services; moreover, from the three measurements, it would be theoretically possible to suppress the anthropogenous noise which occurs on local-scale—the only correlated information is the Earth’s field variation, which is homogeneous enough across the three stations span (100–200 km). From Figure 1, it is evident that mainly the DC railways (which are far away from BDV observatory) will have a detrimental effect on anthropogenic noise on both variometers running at PLM and KEL. The site of Průhonice (PRU) just at the outskirts of Prague is also shown; it served as a geomagnetic observatory station from 1946 to 1967, when the magnetic observations have been moved to a much quieter location of Budkov.
2. Instrument Setup and Site Limitations of PLM and KEL Stations
The station of Polom (PLM) has been in service since late 2016. The site is a property of Czech Army and is being run in collaboration with the Institute of Geophysics, CAS and provides important seismic, meteorological, and geodetic data [8]. CTU and IG CAS took the opportunity to install a fluxgate variometer instrument [4] in the already magnetically prescreened and prepared locality. Because of the recent installation, the temperature-stabilized hut is not yet available, and therefore, the variometer sensor and also electronics are operating at ambient temperatures, although protected from the elements. The site is also equipped with nonmagnetic pillars for obtaining “absolute” magnetic measurements, i.e., measurements of local inclination and declination hand-to-hand with total field intensity, which are usually obtained using a portable Overhauser magnetometer [9].
The station of Kelčany (KEL) is privately owned and is being run by the members of Magnetic Laboratory at the Department of Measurement, FEE CTU Prague. The advantage of the site compared to PLM is the underground location of the sensor and electronics (in a dual-purpose wine-cellar), which allows for less than ±5°C yearly temperature variation. Careful magnetic mapping has been done before installation, the site was cleaned of ferromagnetic objects and a nonmagnetic pillar for the instrument was built. The site is running since 2015 and is, advantageously, on roughly the same latitude as the BDV observatory. The data are publicly available [10].
The magnetometers installed both at KEL and PLM stations were manufactured at the CTU using low-noise race-track fluxgate sensors, exhibiting ~20 pT digital noise floor and showing high geometrical and temperature stability. The triaxial sensor head at PLM is moreover made from MACOR machinable ceramics; the head is further fixed on a marble plate [3]—see Figure 2. Neither KEL nor PLM facilitate a scalar magnetometer; thus, the total field data are calculated from the three orthogonal field components.
(a)
(b)
After initial trials at the PLM station, where the ambient temperature can change from −20°C to +40°C, we decided to orient the sensor to the “UVZ” orientation [11]. This means that the two horizontal axes are oriented ±45° from local meridian—in this manner, both horizontal axes are measuring roughly the same magnetic field (about 15,000 nT at our location). The NEZ or HDZ components are computed numerically [9], so the offset drifts and mechanical instability in azimuth are of less significance than if measuring the E or D component directly. Also, the compensating current in all axes is large enough (few mA) not to be influenced by cable leakage currents. The UVZ orientation is also beneficial for obtaining a simple thermal drift model as shown later.
A comparison of anthropogenous noise observed at BDV, KEL, and PLM stations is shown in time-domain in Figure 3, where the calculated total field from both KEL and PLM vectorial readings is compared to total field measurements at BDV observatory provided by Overhauser magnetometer.
(a)
(b)
The anthropogenous noise at both PLM and KEL is larger than at the BDV observatory (about 0.2 nT peak-peak); however, at KEL also isolated peaks occur with an amplitude up to several tens of nT p-p, which have been later traced down to be caused by nearby car traffic, as further shown in Section 2.1.
Figure 4 further shows aligned spectrograms from three days of 1-second data obtained at all three sites (21–23/7/2017). At PLM, the clean nights are alternated with noisy daytime periods due to the ~40 km distant DC railway and light urban rail. Although the noise at BDV station is very low, the used instrument (DMI fluxgate variometer) has large intrinsic noise, so actually, during the quiet night periods (with almost no electric train traffic), the PLM data are less noisy due to the used variometer. The KEL data on the other hand suffer from increased anthropogenic noise even during the night, since the sensor is located in a residential location—the broadband daytime noise is another 10 dB above PLM. Another 10–20 dB noise increase in short bursts has been traced down as local car traffic.
2.1. Car Traffic and Magnetic Noise
To confirm the origin of the excess-noise at KEL site, the passing cars (the local street is about 25 m away from the sensor location) have been observed by a web-camera and by a cellphone video recording, respectively, and compared to magnetic data—a sketch displaying the actual setup at KEL site is shown in Figure 5.
To be able to detect, mark, and possibly remove the passing car’s magnetic signature, an axial (dBy/dy) fluxgate gradiometer has been created in N-S direction by placing a second sensor coaxial to the variometer head. This second sensor has been placed approx. 5 m away from the variometer (position 1), closer to the street; later also a short-baseline gradiometer was placed at position 2. The peaks obtained from the axial gradient data correspond with the peaks of the variometer data; however, the gradiometer noise floor is still high for detecting spikes less than about 5 nT p-p.
In Figure 6, the vehicle occurrences have been drawn into the magnetic field recorded. Axial gradient (in N-S direction) and magnetic field (N-S component) are shown—the recorded spikes are in the order of tens nT p-p (even larger for vans/busses); we show that there is a clear correlation of the spikes and car traffic.
3. Methods of Improving the Real-World Performance at KEL and PLM Stations
At KEL and PLM sites, we are experiencing two difficulties: temperature drifts due to seasonal and diurnal changes of the ambient temperature (PLM) and large noise due to car occurrences (KEL). It should be noted that the car-induced spikes cannot be simply filtered out with a low-pass filter since the peaks in the individual axes are not bipolar, thus any low-pass filtration would introduce artefacts in the measured data. Moreover, for a strong magnetic source (e.g., a bus-coach), the disturbance occurs even for 20 seconds.
3.1. Correcting for Temperature Drifts (PLM)
It is obvious that if the temperature coefficients of the sensor are known (i.e., offset and gain temperature coefficient), one could recalculate and obtain drift-free data. However, the temperature drift in a fluxgate magnetometer (assuming that the electronic is drift-free) is caused by multiple effects [12], e.g., by the temperature of the excitation tank capacitor, by the dimensional expansion of the feedback/pick-up coil, or due to the expansion of the triaxial holder material and its base. Moreover, it is difficult to calibrate the whole setup as the sensor, and its base are quite bulky. As we have selected the UVZ orientation of the sensor, the ~20–30 ppm⋅K−1 sensitivity drifts dominate in all axes, simplifying further modeling of thermal response (With NEZ orientation, the E axis drift would be dominated by the offset drift which is a combination of electronic and sensor drifts; however, the electronics and sensor head are in our case at different positions and temperatures and exhibit different thermal mass. With HDZ orientation, the D component (~2000 nT) would be influenced both by the offset and gain drifts. In both NEZ and HDZ cases, also mechanical directional instability would have to be modeled.). The predicted sensitivity drift is 0.4 nT⋅K−1 for each horizontal axis and 1.1 nT⋅K−1 for the vertical axis, respectively. Utilizing a “full-field” variometer, thus measuring in a feedback loop all the three vector components of the magnetic field at once, allows us for calculation of the total magnetic field (the scalar vector magnitude). Both variometers at PLM and KEL have been calibrated with the “scalar method” [4] for their offsets, gains, and orthogonalities, so the only difference to a drift-free scalar measurement from an Overhauser magnetometer are then the magnetometer drifts itself.
We did this for the PLM variometer by comparing the Overhauser readings obtained at BDV observatory to the calculated total field from PLM—in this case, we assume that on the local scale, the measurements at the two localities, which do not exhibit geologic anomalies, will differ only by a stable offset BOff. This was also verified during multiple onsite measurements with an Overhauser magnetometer at different times and temperatures (we could not yet perform a long-term scalar measurement due to the lacking infrastructure).
To find the actual variometer drifts in all three axes, we have utilized a least-squares fitting method, which generally minimizes the difference BDiff between the scalar reading BBDV at the BDV observatory and the calculated scalar value at PLM from the three individual components B1, B2, B3. Thus we try to minimize BDiff from a large set of following equations:
The solution of equation (1) was found with a constrained fminsearch function in MATLAB R2015 [13], and the offset agreed well with the one obtained from onsite Overhauser measurements. After correcting on the obtained drift constants α, β, and γ [T⋅K−1], we were able to largely suppress the temperature drifts in all three axes. The dataset we have used was from February 2018, which allowed for large temperature span between +17°C and −12°C, see Figure 7.
(a)
(b)
We could improve the results even further by introducing a lag of 800 s which was experimentally obtained by calculating the cross-correlation between the total field differences and temperature—this delay is believed to originate from the fact that the temperature measurements occur at the MACOR cube where the sensors are located, but significant part of the drifts can be caused by the excitation capacitor temperature coefficient [14]—the capacitor is heated only by radiation, since it is thermally connected to the MACOR cube only by its thin leads. After introducing this delay, the calculated values were the following:
The value of α roughly corresponds to the 30 ppm predicted drift (vertical axis measuring approx. 44,000 nT). Also, the γ value corresponds to an expected value for a horizontal sensor. However, the obtained value of β is unexpected, since both horizontal sensors should exhibit the same values or at least the same order of magnitude. We cannot currently offer other explanation than a faulty sensor deployed at this position.
After rotating the temperature-compensated PLM vector readings with a 3 × 3 matrix, which reorients the sensor at PLM to the orientation at BDV, we were able to show that the temperature compensation was successful also in the individual components, see Figure 8.
3.2. Suppressing the Car-Induced Disturbances (KEL)
As we have shown in paragraph 2.1, there is a clear correlation between the disturbances occurring at KEL and the car traffic. Thus we decided to create a simple model, assuming the following simplifications: (1)The car at the ~25 m+ distance can be well modeled as a single magnetic dipole(2)We neglect the road curvature and assume it in E-W direction (x-axis)(3)The magnetic moment magnitude and orientation are stable during the car passage, since it keeps its orientation to the Earth’s magnetic field(4)The occurrence of the maximum axial gradient in the y-axis, which occurs when the car radial distance to the sensor is smallest, defines the symmetry point of the car movement(5)The maximum axial gradient in the y-axis (N-S component) occurs defines symmetry of the car movement(6)The car does not change its speed significantly(7)In 10 seconds, the car is distant enough not to give any significant (>0.1 nT) disturbance(8)During the ~20 s car passage, the Earth’s magnetic field changes only linearly(9)The occurrences of cars in both lanes are not frequent.
Thus when fulfilling the above assumptions, we can write for the vectors of observed field and vector of car disturbance field utilizing the well-known equation for magnetic field of a magnetic dipole with magnetic moment position vector and an (orthogonal) rotation matrix : where the position vector size (radial distance) is calculated as
The position vector coordinates of equations (3) and (4) are aligned with the magnetic moment coordinates of the dipolar source (hence the need for rotational matrix to align with the coordinates). However, we are not interested in the real orientation of the magnetic moment vector in this case, so there is no need to calculate for , which would further complicate the problem.
To find the “true, disturbance-free” of equation (3), we have implemented a least-squares fitting-based algorithm, which relies on reading from a gradiometer placed close to the street, allowing for detection of the passing cars through the “point of symmetry”. In other words, we are trying to fit the magnetic field at every sampling point using equation (3). To achieve this, we implement a set of equations which describe the magnetic field during the short disturbance, which occurs due to the moving magnetic moment. The dataset for the optimization is centered at the gradient peak and is usually 10 + 10 seconds long (Figure 4). Since we assume a simple trajectory as of Figure 4, we are trying to find position components x and y, whereas the only changing is the x, since y is constant and equals zero in our case. Due to the constant sampling time, we can express x as a linear function of time and vehicle speed , both of which we assume constant. Since the sensor at KEL is sampling at 206.5 samples⋅s−1, there are enough equations during the car passage, even after FIR filtering of the data (to remove 50 Hz mains disturbances) and smoothing. The optimization result of equation (3) is then the “true” Earth’s field vector , the magnetic moment vector , the car speed , the initial position x0, and the constant y distance together with the time-derivative of the Earth’s magnetic field during the fitting period.
The optimization started only when the gradient amplitude in the observed interval crosses a preset threshold in order to run only for disturbances significantly larger than overall system noise. In our case, the threshold has been set to 2 nT⋅m−1. The algorithm also contains bounds and tests to compensate only using the expected values (car speed ~ 5–20 m⋅s−1, y between 20 and 30 meters, fitted |m| below 600 A⋅m2). For details of the algorithm and used functions, see Appendix.
3.2.1. Webcam-Trial: Magnetic Moment Statistics
We applied the compensating algorithm on the dataset from the verification video trial (see Chapter 2.1) in order to check the feasibility of our model. In the 50-minute dataset, we had ~270 car passages: 2 bus-coaches, 4 vans, and 17 motorbikes and the rest were passenger cars (see Figure 6 for the induced disturbances); they were evenly distributed in the close and far lane (130 vs. 133 occurrences). The resulting “typical” magnetic moment for passenger cars and busses was found as 250 ± 50 A⋅m2 and 520 ± 50 A⋅m2, respectively (the compensating algorithm did not start for the motorbikes due to the gradient threshold). In Figure 9, we show statistics of the individual magnetic moment components (only passenger cars shown), from which it is evident that the largest component is the vertical one which tends to be oriented along the magnetic flux lines, i.e., the component does not change its sign depending on car trajectory orientation.
(a)
(b)
(c)
3.2.2. Compensating the Disturbances
An example result for compensating a single car disturbance is shown in Figure 10—the original data, the fitted dipole from moving the car, and the data after disturbance compensation are shown for x, y, and z field components. The noise clearly decreased in the cleaned-up dataset—see Figure 11 for a close-up of time-domain and Figure 12 for two spectrograms of 14-hour data.
(a)
(b)
(c)
(a)
(b)
(a)
(b)
4. Conclusion
We show an approach on how to deal with (1) unstabilized ambient temperature and (2) man-made disturbances occurring at two variometric stations. Whereas the first problem is usually solved by temperature stabilization—active or passive [15]; our approach shows that a “moderate” performance can be achieved even at ambient temperatures. Although the presented method seems straightforward, we are not aware of other ambient-temperature magnetic field stations utilizing such long-term calibration and compensation. The overall maximum residual drift of 5 nT p-p during 30°C temperature swing was achieved, which is even in accordance with INTERMAGNET standards [16], where the instrument should keep 0.25 nT⋅C−1 for vectorial readings, but it still does not fulfill the required 1 nT accuracy for scalar values. However, for our purposes, this approach brings fast and reliable results as we can really choose the ambient-run site at PLM as a redundant source of magnetic data. The steps and results shown here can be beneficial to many “repeat” stations, which usually run at ambient temperatures and which are supplementing the magnetic observatories. Even better results can then be expected if the sensor is, i.e., buried at 1-2 meters to avoid such large temperature fluctuations, and of course, when utilizing at least a moderate temperature stabilization (±2°C), the residual drifts after fitting would be one order of magnitude less than those presented.
As for the second problem, fitting and cleaning of a 14-hour 1-second dataset took less than 60 seconds on a Core-i7 PC using MATLAB R2015, so off-line postprocessing of daily data could be viable even in embedded systems running Linux and using Python fitting libraries. The fitting speed and accuracy can be improved by having an a priori knowledge of the target y-distance (i.e., knowing the lane) and/or the car speed. For achieving this, distance calculation from axial gradiometer reading [17] or speed measurement with magnetic sensor [18] could be used. Also, improved detection of the passing cars can be facilitated with a short-baseline gradiometer placed as close to the street as possible, which also results in a better approximation of the gradient by the calculated field difference. As for the problem with two cars passing in adjacent lanes, it could be possibly solved by employing a second order model with two car trajectories and vector summation of the magnetic signature; for detection of this situation could symmetrically deploy two short-baselines along the street and observe their output in time.
We again emphasize that this method, when compared to the usual ways—disregarding or interpolating the data heavy low-pass filtering [9, 19]—allows not only to keep the true DC value but also does not destroy the “high frequency” component originating from various physical phenomena (field oscillations, solar storm onsets, etc.).
Appendix
The basic MATLAB function utilized in the optimization problem for car-induced disturbances cancelling is shown in this appendix. The field-fit function enters the least-squares curve-fitting algorithm lsqcurvefit of the Optimization Toolbox (MATLAB 2015) with a guess of the values and lower and upper bounds, where xdata is the time vector and BObs is the vector of the observed magnetic field. The function dip2B to calculate the field from a point-like dipolar source follows the dipolar term in equation (3); an example implementation can be accessed, e.g., in [20].
| ||||||||||||||||||||||||||||||||||||||||
Data Availability
The geomagnetic data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare no conflict of interest.
Authors’ Contributions
Conceptualization, investigation, software, methodology, data curation, and visualization was provided by M. J.; writing was performed by M. J., M. V., and T.B.; resources and data collection were provided by M. B., M. J., M. V., and T. B.
Acknowledgments
We would like to thank the Military Geographic and Hydrometeorologic Office of the Czech Army for the possibility to work and acquire data at the Meteorologic and Seismic station Polom, including the complimentary accommodation. We would also like to thank Jaroslav Pokorný for his logistic and technical support during onsite measurements.