Network of sensitive magnetometers for urban studies

The magnetic signature of an urban environment is investigated using a geographically distributed network of fluxgate magnetometers deployed in and around Berkeley, California. The system hardware and software are described and results from initial operation of the network are reported. The sensors sample the vector magnetic field with a 4 kHz resolution and are sensitive to fluctuations below 0.1 $\textrm{nT}/\sqrt{\textrm{Hz}}$. Data from separate stations are synchronized to around $\pm100$ $\mu{s}$ using GPS and computer system clocks. Data from all sensors are automatically uploaded to a central server. Anomalous events, such as lightning strikes, have been observed. A wavelet analysis is used to study observations over a wide range of temporal scales up to daily variations that show strong differences between weekend and weekdays. The Bay Area Rapid Transit (BART) is identified as the most dominant signal from these observations and a superposed epoch analysis is used to study and extract the BART signal. Initial results of the correlation between sensors are also presented.


I. INTRODUCTION
The study of fluctuating magnetic fields has found widespread application in a variety of disciplines. High-cadence (or effective sampling rate) magnetic field measurements have been used for magnetic-anomaly detection (MAD), which has numerous security applications, such as naval defense and unexploded ordnance detection 1 . Additionally, geographically distributed magnetometers, operating at high sample rates, have been designed to study the global distribution of lightning strikes using the Schumann resonance 2 .
Low-frequency magnetic field measurements, typically corresponding to large length scales, provide important information relating to the nature of magnetic sources in the Earth's core, maps of near-surface fields, as well as crustal composition and structure models 3 . An international consortium of magnetometer arrays, known as SuperMag, comprises approximately 300 magnetometers operated by numerous organizations 4,5 . The magnetic field data collected from these stations is uniformly processed and provided to users in a common coordinate system and timebase 5 . Such data are important to global-positioning-system (GPS)-free navigation, radiation-hazard prediction 6 , climate and a) Electronic mail: tbowen@berkeley.edu weather modeling, and fundamental geophysics research. Additionally, measurements of the auroral magnetic field are necessary in testing models for space-weather prediction, which aims to mitigate hazards to satellite communications and GPS from solar storms [7][8][9] .
Magnetometry has additionally been applied in the search for earthquake precursors. Anomalous enhancements in ultralow frequency (ULF) magnetic fields were reported leading up to the October 17, 1989 Loma Prieta earthquake 10 . Similar anomalous geomagnetic behavior was observed for the month preceding the 1999 Chi-Chi earthquake in Taiwan 11 . Recently, attempts have been made to study precursors to the 2011 Tohoku earthquake in Japan 12 .
Despite its numerous applications, magnetometry is frequently limited by contamination from unrelated, and often unknown, sources. In their search for earthquake precursors, Fraser-Smith et al. (1978) found that magnetic noise from the Bay Area Rapid Transit (BART) system dominated their ULF sensors 13 . This noise is evident in Fig. 1, which depicts the magnitude of the magnetic field recorded at the University of California (UC) Berkeley Botanical Garden. Despite the obvious presence of local permanent or time-varying magnetic contamination, the recorded magnetic field is dominated by fluctuations beyond the expected typical geomagnetic values. These fluctuations diminish dramatically between approximately 1 AM and 5 AM local time, indicating that these are the same fluctuations that are attributed to the operation of BART 13 .
Even during the magnetically quieter nighttime period, there remain fluctuations which exceed expected geophysical values. Certainly, some of the field variation can be attributed to variations in the ionospheric dynamo and other natural sources: similar trends appear in the magnetic record from the Botanical Garden as well as Intermagnet data from a magnetometer located in the nearby city of Fresno 15,16 . Nevertheless, the contributions of human activity to these nighttime fluctuations remain poorly understood.
Only recently have investigators been able to accumulate and study the increasing quantity of spatially and temporally granular data that characterize the evolving state of a city. These data include information from social networks (e.g., twitter), financial transactions, transportation systems, environmental markers (pollution, temperature, etc.) and a wide range of other physical quantities. Passive observations of cities not only serve as a means of quantifying urban functioning for the purpose of characterizing cities as complex systems of study, but they can also yield tremendous benefit to city agencies and policy makers.
Recent work has shown that broadband visible observations at night can identify patterns of light that can be used to measure total energy consumption and public health effects of light pollution on circadian rhythms 17 . High frequency (∼120 Hz) measurements of urban lighting allow for phase change detection of fluorescent and incandescent flicker, which may be used as a predictor for power outages 18 . In addition, infrared hyperspectral observations can determine the molecular content of pollution plumes produced by building energy consumption providing a powerful method for environmental monitoring of cities 19 .
Here we report on the development of a synchronized magnetometer array in Berkeley, California for the purpose of studying the signature of urban magnetic field over a range of spatiotemporal scales. The array, currently consisting of four magnetometers operating at a 3960 Hz sample rate, will make sustained and continued measurements of the urban field over years, observing the city's dynamic magnetic signature. Through systematically observing magnetic signatures of a city, we hope to complement advances made elsewhere in urban informatics and applied physics to provide a deeper understanding of urban magnetic noise for researchers in geophysics and earth science.
In Sec. II we briefly describe the components and performance of the hardware and software implemented in our magnetometer array. We emphasize our preference for commercially available hardware and advanced timing algorithms and utilize techniques similar to those of the Global Network of Optical Magnetometers for Exotic physics searches (GNOME) project 20 .
The analysis of our data is presented in Sec. III. In Sec. III A we present the signatures of several common urban sources and observe the signature of a lightning strike. Station data are analyzed and compared in Sec. III B. Clear variations between weekday and weekend, day and night, and distance from BART tracks are observed. Measurements are examined with wavelets in Sec. III C. Initial results of correlating station data are presented. In Sec. IV, we present an initial method to isolate the BART signal. Conclusions and directions for future research are presented in Sec. V.

II. INSTRUMENTATION, HARDWARE, AND DATA ACQUISITION
A. Description of the system Our magnetometer network consists of several spatially separated magnetometers with high-precision timing for correlation analysis. Figure 2 shows a map of Berkeley with the sensor locations of the network. Each station consists of a commercially available fluxgate magnetometer, a general-purpose laptop computer, and an inexpensive GPS receiver. Our approach avoids bulky and expensive hardware, favoring consumer components wherever possible. As a result, we reduce the cost of the acquisition system and enable portability through battery operation. However, achieving the desired timing precision (∼100 µs) with affordable commercial hardware requires implementing a customized timing algorithm.
A schematic of a single magnetometer is presented in Fig. 3. Each setup is controlled by a computer (PC, ASUS X200M) running the Windows 10 operating system (OS). The PC acquires magnetic field data from the DAQ (Biomed eMains Universal Serial Bus, USB, 24 bit) and timing data from the GPS receiver (Garmin 18x LVC). The DAQ continuously samples the fluxgate magnetometer (MAG, Biomed eFM3A) at a rate of 3960 sample/s. Absolute timing data are provided once per second by the GPS receiver, which is connected through a powered high-speed RS232-to-USB converter (SIO-U232-59). The GPS pulse-per-second signal, with 1 µs accuracy, is routed to the computer through the carrier detect (CD) pin of the RS232 converter. Data from the DAQ arrive in packets of 138 vector samples approximately every 35 ms. As the data are received, they are recorded together with the GPS information and the computer system clock. Data are uploaded via wireless internet to a shared Google Drive folder. FIG. 1. The Earth's magnetic field, as recorded at the UC Berkeley Botanical Garden. The red trace represents measurements taken by an all-optical self-oscillating magnetometer designed and constructed at UC Berkeley in partnership with Southwest Sciences, Inc. 14 . The green and blue traces are data taken with a commercially available cesium magnetometer from Geometrics, Inc. The gray trace represents the average field seen by a geomagnetic observatory located in Fresno, CA.

B. Time synchronization and filtering
Time intervals between the GPS updates are measured by the computer system clock (performance counter), which runs at 2533200 Hz. A linear fit model is used to determine the absolute system time relative to the GPS. Only the GPS timing data from the last 120 seconds are used to determine linear fit parameters. When a magnetic-field data packet is received, the acquisition system immediately records the performance counter value. The packet time-tag is determined from interpolating the linear fit GPS time to the performance counter value. Typical jitter of the inferred time stamps is 120 µs and is limited by the USB latency 21 .
Some of the data packets cannot be processed immediately due to OS limitations. These packets are delayed by up to several milliseconds before being delivered to the data-acquisition software. The number of the delayed packets depends on the system load. During normal operation of a magnetometer, the average fraction of the delayed packets is about 3%. The intervals between both GPS and packet arrival times are measured with the performance counter in order to identify the delayed packets. Any GPS data that arrive more than 50 µs late are discarded from the linear-fit model. When magnetometer data packet arrives with a 200 µs or greater discrepancy from the expected time, their time stamp is replaced with the expected arrival time, which is inferred from the linear fit. Our time filtering algorithm and data acquisition software are publicly available on GitHub 22 .

C. Performance characterization
In order to characterize the performance of the data acquisition system, we apply a reference signal simultaneously to the sensors. All four sensors are placed together into a single Helmholtz coil system driven by a pulse generator. The amplitude of the pulses is 2 µT, the period is 200 ms, and the duty cycle is 50% (Fig. 4a). The top and bottom rows in Fig. 4 represent the data before and after application of the timing-correction algorithm. The inset in Fig. 4a demonstrates how a delay in retrieving a data packet disrupts the timing of the field pulses. When a data packet is delayed, the magnetic-field samples are distributed over a slightly larger time period, which affects the estimated time of the field change. In Fig. 4, interpolated time of the falling edge from sensor four has been several milliseconds after the field changed, causing the zero crossing of the square wave to shift in time from the other sensors. Figure 4b shows the time discrepancies between the interpolated and expected square wave zero crossings. The zero crossings are recorded with 120 µs standard deviation from the expected interval; however, there are a large number of outliers with up to a 10 ms discrepancy. Figure 4c shows the histogram of the data from 4b, where the red curve represents the best fit Gaussian. After the timing correction algorithm is applied to the data, time stamps associated with outlier packets are replaced with the expected arrival times (Fig. 4 d,e,f). The remaining jitter is a Gaussian distributed with an error of 120 µs. Power spectral densities were calculated for each interval, with the noise floor taken as the ensemble average of the interval spectra. The noise floor varies between individual axes, with the most noise observed on the Z-axis. For all three directions, the noise floor is constant between ≈ 2 Hz and ≈ 700 Hz. Narrowband spectral features from 60 Hz and harmonics are easily observed in the data. For frequencies above 1 Hz, the noise floor is uniformly below 0.1 nT/ √ Hz. The peak observed in the noise floor, predominately in the Z and Y channels, between ∼ 300 and ∼1500 Hz, is likely due to the operation of the fluxgate electronics.

A. Observations of urban magnetic signatures
The portability of our sensors has enabled the direct measurement of several urban field sources. Figure 6 shows the magnetic signatures of several field sources associated with transportation: traffic on a freeway, as well as both Amtrak and BART trains. Figure 7 shows a spike in the magnetic field due to a lightning strike recorded by three geographically distributed sensors. This lightning strike, observed before the implementation of the timing algo-rithm, highlights the need for time corrected data: e.g. interpreting the time lag between sensors as a physical transit time between stations corresponds to unreasonable spatial separations of ∼15, 000 km. Unfortunately, no further lightning strikes have been captured since im-  13 . The dependence of these fluctuations on proximity to the BART lines, and their correspondence with train timetable, led the authors to attribute this ULF signal to currents in the BART rail system. Subsequent ULF measurements made in 23 at a distance of 100 m of BART suggested periodic bursts of magnetic field at roughly the periodicity of the BART train. The location of our network stations (up to 2 kms from BART line) and length of recorded intervals ensure that urban signatures, such as BART, will be present in our data.
The timing algorithm provides sub-millisecond resolution for MAD; however many magnetic signatures related to anthropogenic activity occur at low frequencies where GPS alone provides adequate timing. To investigate urban magnetic fluctuations, we decimate the full 3,960 sample/s to a 1 sample/s cadence. Antialiasing is accomplished with a moving average (boxcar) filter. Low-cadence data provide adequate resolution for correlating multi-station observations, while simulta-neously removing the 60 Hz power signal. We only use data from Stations 2, 3 and 4 in multi-station comparisons. Station 1 served mainly as an engineering unit for characterisation measurements and other testing. A magnetically quiet period corresponding to BART non-operating hours 24 is evident in the records of the three active sensors (stations 2-4). Figure 9 shows the distribution of magnitude fluctuations binned at 10 −3 nT. The record naturally separates into two intervals, corresponding to the hours when BART trains are running (roughly, from 7:55AM to 1:26AM) and hours when BART trains are inactive. We refer to these intervals as "daytime" and "nighttime", respectively. The magnetic measurements reveal characteristically different distribution functions in daytime and nighttime. For  The two vertical bars in Fig. 10 mark the period for which BART is inactive. Increased field variance, i.e. signal power, is clearly tied to the operation of BART. Figure 10 also highlights the relatively constant daytime variance in each station. Additionally, there appear to be approximately constant ratios between the daytime variance (power) observed by each sensor. This suggests that the background noise level observed in each station is set by the distance to BART. Identifying other anthropogenic fields (for example, traffic) is complicated by the large BART signal. Accordingly, identification of further urban signals requires a thorough characterization of the magnetic background generated by BART.

C. Time-Frequency (Wavelet) Analysis
Frequency-domain analysis is typically used to reveal the spectral composition of a magnetic time series. Localizing the distribution of spectral power in time requires simultaneous analysis in both time and frequency domains. We implement a continuous wavelet transform (CWT), using Morlet wavelets 25 . to investigate the time-frequency dis-tribution of low-frequency fluctuations associated with BART. The unnormalized Morlet wavelet function ψ(τ ) is a Gaussian modulated complex exponential, with non-dimensional time and frequency parameters τ and ω 0 . A value of ω 0 = 6 meets the admissibility conditions prescribed in Ref. 26 and is commonly used across disciplines 25,27 . At each time step, the CWT of the magnetic field B is defined by the convolution of the time series record with a set of scaled wavelets, which are normalized to maintain unit energy at each scale. The CWT provides a scale independent analysis of time-localized signals, and is additionally insensitive to time series with variable averages (non-stationary signals). These qualities provide some advantage over alternative time-frequency analysis techniques, such as    the windowed Fourier transform, which calculates the spectral power density in a sliding window applied to the time series. Introducing a window imposes a preferred scale which can complicate analysis of a signal's spectral composition. For example, low-frequency components, with periods longer than the sliding window scale, are aliased into the range of frequencies allowed by the window, thereby degrading the estimate of spectral density.
Full day, 1 sample/s cadence, wavelet power spectral densities |W (s, t)| 2 (µT 2 /Hz) for stations 3 and 4 on Sunday, 03/20/2016 (PST) are displayed in Fig. 11. These spectrograms prominently display the quiet night-time period. Additionally, strong power is observed in several scales corresponding to a fundamental 20-minute period (8.33×10 −4 Hz) and associated higher harmonics. This 20-minute period coincides with the Sunday BART timetable on the geographically closest BART line (Richmond-Fremont). The black lines display the region where boundary effects are likely -this region, known as the cone of influence (COI), corresponds to the e-folding time for the wavelet response to an impulse function. This plot immediately suggests that stations 3 and 4 are highly correlated, while station 2 is anti-correlated with stations 3 and 4. Indeed, this is verified by the bottom panel of Fig. 12, which shows the cross correlation coefficients C ij (τ ) calculated for the 24 hour time series: where N is the record length, τ is a translation between time series, and n is the sample index 28 . It is clear that stations 3 and 4 are in phase, while station 2 is out of phase of the other two instruments. These phase relationships which correspond to the geographical location of the sensors on the east/west side of the BART line (stations 3 and 4 are located east of the rails, while station 2 is located to the west, c.f. Fig. 2), suggest that the magnetic field generated by the BART has a strong azimuthal symmetry around the BART rail. Future work will look to determine the multipole components (e.g. line current, dipole, and higher corresponding to the BART field. Figure 13 shows full day wavelet spectral densities for station 2 on both Wednesday, 03/16/2016 and Sunday, 03/20/2016. The quiet BART night is much shorter on Wednesday; this corresponds with the different weekend and weekday timetables. Additionally, Fig. 13 demonstrates the absence of a strong 20-minute period in the Wednesday data, instead revealing more complex spectral signatures. The increase in complexity observed in the weekday spectrogram is most certainly due to increases in train frequency, the addition of another active BART train, and variability associated with commuter hours. In our future work we will fully explore the effect of train variability on the urban magnetic field using correlated observations from the entire network. The measurements of 23 made between 0.001 and 4 Hz, at a range of several hundred meters from the BART rail, captured bursts of magnetic field corresponding approximately to the train schedule, but with an irregular variation in the waveform of the magnetic field. In contrast, our observations reveal a highly regular signature, with multiple spectral components, occurring at the BART  30 . We identified 46 sharp peaks in the magnetic field occurring with an approximately 20minute period (e.g. Fig. 8). From these 46 peaks, an ensemble [X(t) i ] of intervals is constructed comprising the 3 minutes preceding and 17 minutes succeeding each individual peak. Averaging over the ensemble of intervals X(t) = i X i (t) reveals a coherent signature with an approximate 20-minute period, Fig. 14(a). The periodic signal observed in the data has the form of a sharp discrete peak of ≈1 nT, followed by an oscillation with a period on order of several minutes. A quantitative comparison is obtained through computing the Pearson correlation of the extracted signal with each interval in the ensemble. On average, the correlation between the extracted signal and observed data has a correlation ofρ = 0.7, with ρ i ranging from 0.1 to 0.85. We can interpret these values as the fraction of power in each interval derived from the  average signature. Figure 14(b) demonstrates high correlation between the extracted average signal with an hour of observations taken from 9-10 AM (PDT). Through extracting periodic magnetic signatures of BART, we enable the identification of transient events associated with BART operation as well as other urban phenomena. Fig.14(c) shows the occurrence of a local magnetic anomaly occurring at 12PM (PDT) in Station 2; in this case ρ = 0.67. The traces from the other stations, suggest that the event observed in station 2 is a local anomaly, not associated with BART. Additionally, Figure 14(d) demonstrates an interval of data (with ρ = 0.17) which includes a global transient feature, likely due to some variation in the BART system. Measurements from a single sensor allow us to identify events which deviate from the correlated periodic observations; our future work will employ the full network of magnetometers to identify correlated signals in both space and time, allowing for an extraction of the magnetic field local to each sensor from the global field dominated by the BART signal.

V. DISCUSSION
An array of four magnetometers has been developed with bandwidth of DC-kHz and sensitivity better than 0.1 nT/ √ Hz. The array is currently deployed in the area surrounding Berkeley, CA, providing measurements of an urban magnetic field. This array is sensitive to both natural magnetic activity, such as lightning and the low frequency variations in the Earth's geomagnetic field, as well as a variety of anthropogenic sources: currents associated with BART, traffic and 60 Hz powerlines. The operation of BART dominates the urban magnetic field generated broadband noise. In addition to this broadband noise, the network has identified the presence of coherent narrowband spectral features originating from the BART. Significant variation in the spectral features is observed between weekends and weekdays corresponding to variations in the BART train schedule. During the hours in which BART is non-operational, the anthropogenically generated fields are significantly decreased and agreement with the USGS magnetic field measurements is observed. However, the nighttime field still contains a number of features not attributable to geophysical activity. Further study is required to determine the nature and sources of these features.
Cross-correlating the sensors at high frequencies requires a high-precision timing algorithm to combine the absolute time, acquired through GPS, with the high-precision computer performance clock local to each station. This algorithm additionally corrects for latency issues associated with the USB interface between the data-acquisition hardware and the laptop operating system. This timing algorithm has been tested using magnetic fields generated by Helmholtz coils. We intend to use the impulsive globally-observable fields generated by lightning to further test our timing algorithm. Our high precision timing will allow for such magnetic anomaly detection on the order of ≈ 100 µs.
This paper presented a proof-of-concept deployment of what, to our knowledge, is the first synchronized network of magnetometers specifically designed for observing the effects of human activity on the magnetic field in an urban environment. Numerous potential applications and directions for future work have emerged. Further development of algorithms to remove the BART signal (or any other dominant signal whose source has been identified) must be developed. These algorithms may need to take advantage of other data sources (e.g., the realtime BART schedule) and machine learning techniques. The study of high frequency response (60 Hz and above) has not yet been pursued. We note that anthropogenic fields mask geophysical field fluctuations, and that study of the latter is facilitated by understanding of anthropogenic noise. The magnetic fields due to humans may reflect identifiable aspects of urban dynamics (beyond BART) and these may have correlations to other measures of urban life (energy consumption being one of the first to consider). Studies of magnetic field correlations and anomalies may be used to identify and study local phenomena (traffic, elevators, etc.). One spin-off of this research may be improved identification and reduction of anthropogenic noise in geomagnetic measurements located in or near urban environments. The ultimate utility of the magnetometer array as an observational platform for urban systems will only become clear with further studies.