Airborne Gravity Data Denoising Based on Empirical Mode Decomposition: A Case Study for SGA-WZ Greenland Test Data

: Surveying the Earth’s gravity field refers to an important domain of Geodesy, involving deep connections with Earth Sciences and Geo-information. Airborne gravimetry is an effective tool for collecting gravity data with mGal accuracy and a spatial resolution of several kilometers. The main obstacle of airborne gravimetry is extracting gravity disturbance from the extremely low signal to noise ratio measuring data. In general, the power of noise concentrates on the higher frequency of measuring data, and a low pass filter can be used to eliminate it. However, the noise could distribute in a broad range of frequency while low pass filter cannot deal with it in pass band of the low pass filter. In order to improve the accuracy of the airborne gravimetry, Empirical Mode Decomposition (EMD)


Introduction
The Earth's gravity field, which is responsible of the structure and shape of the earth, is of primary importance for much Earth science based research and applications as well as for practical applications in geo-information.Airborne gravimetry, which is able to determine medium to high resolution gravity, has an advantage of providing a quick coverage of large areas of the Earth with significant economy over other methods [1][2][3][4][5].Over the last two decades, airborne gravimetry had a great development to the point that advanced strapdown systems can map the gravity field at spatial resolution of several kilometers with the accuracy level of a few mGal [6][7][8].
The aim of the merged measurement system in airborne gravimetry is to estimate the gravity disturbance vector [9].According to Newton's second law of motion and gravity disturbance's definition [10], in an inertial frame, the gravity disturbance vector can be reckoned directly from the total acceleration f i sensed by accelerometers of the gravimeter, the kinematic acceleration ̈ of the aircraft and the normal gravity field ϓ.The relationship can be expressed by the following equation [11] Seemingly simple, the practical implementation of Equation ( 1) is quite involved-the noise of the system is highly large.The noise to signal ratios of the measured gravity data without filtering can be 1000 mGal and more.Much of the noise belongs to the high frequency noise caused by aircraft vibration and the expansion of GPS observation noise in the process of computing aircraft acceleration [6,9].The obvious way to eliminate the high frequency noise is a low-pass filter.The finite impulse response (FIR) and infinite impulse response (IIR) filters are the two widely used low-pass filters in airborne gravity data processing.For instance, the University of Calgary has used the FIR filter to process of strapdown INS (SINS) data [7] and Sun has tried the FIR filter at extracting gravity disturbance of the LaCoste and Romberg S-type gravimeter (LCR) test data in China [12,13].The Butterworth filter which is one kind of IIR filter, has been treated by Forsberg for processing LCR test data in Greenland and United Arab Emirates [14,15].Some other low pass smoothing method used in airborne gravimetry can be found in Jekeli [16], Kwon [11], Senobari [17], Studinger [18], etc.The obvious advantage of low-pass filter is simple and easy implementation.However, it is fact that high frequency content does not necessarily belong to noise only, especially for the airborne system, which is in a high dynamic environment.Moreover, noise could be present over a broad range of frequencies.The effects of attitude errors and aircraft dynamics often evoke observation errors at lower frequency; low-pass filters will not work except if the length of the filter is no less than the period of phugoid motion [19,20].Even though both noise and signals are eliminated in the lower wavelength, some noise still exists in the pass-band which will contaminate the accuracy of the result [5].
The empirical mode decomposition (EMD) proposed in 1998 by Huang of NASA has the potential of estimating and analyzing the trend of the data [21].It is a time-frequency analysis method of processing a nonlinear non-stationary signal, which has been applied in the fields of nonlinear system analysis, geophysics, meteorology and biomechanics [22][23][24].However, there are few papers that have applied EMD for denoising raw airborne gravity data.Although Lu [25] and Cai [26] have attempted to denoise airborne gravity data by EMD, more research is needed to test the efficiency of EMD.
In this study, we try to apply the EMD for denoising the raw airborne gravity data obtained by the strapdown airborne gravimeter SGA-WZ in Greenland in 2012.Firstly, the preliminaries on EMD denoising are shown.After descripting and analyzing the data obtained from the primary repeated flight profiles of SGA-WZ, EMD is employed to estimate gravity disturbances from the raw data.

Principle of EMD Denoising
EMD, which has the potential of estimating the trend of data is a new technique of decomposing a given signal into a set of limited Intrinsic Mode Functions (IMFs) [21].The IMF is one kind of function in which the number of zero-crossings equals to the number of extrema, and the upper and lower envelopes defined by the local maxima and minima respectively are symmetric [21,27].
Given any discrete time signal x(n), n = 1, 2, 3, …, N (N is the total sample number), the extracting process of IMFs from the signal x(n) is an iterative procedure called sifting algorithm.Figure 1 gives the main process of the sifting algorithm which consists of the following steps: Step 1: Let r(n) = x(n) as the initial step.
Step 5: Get the residual by subtracting e(n) from the signal: . Check the residual r(n) and whether it meets the requirements of an IMF.If the residual doesn't satisfy the definition of an IMF, repeat steps from step 2 to step 5 until it does; if r(n) is an IMF, then the first IMF is obtained which is imf(n) = r(n).
Step 6: After obtaining a new IMF, remove r(n) from the original signal x(n) and get the new signal under examination is x(n) = x(n) − imf(n).Repeat previous steps and the sifting process will stop when the final residual becomes a monotonic function.
Finally we achieved the signal expressed as the result of EMD: in which K is the total number of IMFs and r(n) is the final residual.Figure 2 shows, for instance, the top six decomposed results of a set of raw data that will be processed in this study.We can find that the noise dominates the higher frequency temporal modes while the signal dominates the lower ones.
According to this phenomenon, it is possible to extract a signal from noisy data.Denote  ̃ as the corrupted IMF with the noise   and   as a pure IMF, so  ̃ can be written as  ̃ =   +   .The goal of denoising is to extract pure IMF from the noisy one.However, in practice, we can only obtain the estimation of the pure IMF  ̂ .Then, the denoised signal can be written as and now the main problem is how to get  ̂ .Since the noise dominates the higher frequency temporal modes while the signal dominates the lower ones, the simplest method is to dispel the higher frequency IMFs.If this happens, the EMD denoising procedure will act like a low pass filter [28].However, the noise can distribute in both high frequency and low frequency IMFs.In this paper, we tried to use a soft thresholding to estimate all IMFs.The thresholding function widely used in wavelet shrinkage denoising is defined as follows [29]: where   is the thresholding of   .Choosing the value of the thresholding is a fundamental problem in order to avoid over-smoothing or under-smoothing.The thresholding value chosen in this study is the universal thresholding which equals √2 ln  [29,30].

Test Description
The data used in this study were from the primary flight profiles of SGA-WZ in central East Greenland in 2012.The purpose of this test, carried out in cooperation between Technical University of Denmark (DTU) Space and National University of Defense Technology (NUDT), is to evaluate repeatability as well as accuracy of the new airborne gravimeter in Arctic regions.
The SGA-WZ together with the GPS receiver was equipped on a Nordlandair Twin Otter in Iceland.To perform the DGPS positioning, a base station was also operated on the ground.Both GPS data were collected at a rate of 1 Hz, and raw accelerometer and gyroscope data from SINS were recorded at a data rate of 2 KHz.
Figure 3 shows the flight profiles in Greenland test and existing ground, marine and airborne coverage of the region.The two primary repeated marine profiles from north to south (shown in bold black lines), which were coincident with the earlier DTU Space flight profiles [14], were in a region with marine data coverage.The total length of each repeated profile was about 300 km, the average flying speed was about 250 km/h and the flying altitude was around 360 m.More details on this campaign could be found in Zhao [31].In order to check the external accuracy of EMD denoising, former airborne gravity data and marine data around this area were used.The former airborne profiles have been flown in 2001 with an LCR gravimeter.The accuracy of the measuring data was no more less 2 mGal under the resolution of around 6 km [14,32].Since the old profiles do not coincide exactly with the new ones, the marine data around this area were used to fill the gaps (See the color lines in Figure 3).First, the marine data was spatially interpolated to the ground track of the flights.Then the marine data numerical combination with the old flight data was therefore used as reference data for the SGA-WZ test.Moreover, the former result from a FIR filter is provided for easy comparison of internal or external accuracy.Note that the resolution and accuracy of the FIR result vary with the length of the filter.For instance, in this case, the resolution of a 120 s (200 s) length FIR is 4 km (7 km), and the internal accuracy is 2.4 mGal (1.5 mGal) for Line A. Here, only the result of a 160 s length FIR is given since the correspondence resolution is same as LCR.

Raw Data Preparation
Equation (1) shows the relationships between the gravity disturbance and the quantities derived from airborne gravimetry system, but it is more convenient to implementing it in practical in a navigation frame as follows: For scalar gravimetry, only the vertical quantity of the gravity (the third component of Equation ( 5)) is considered which can be written as where subscripts stand for Down in a local-level ellipsoidal frame,   is also called the Eötvös correction [7,33].Figure 4 shows a summary of the process of using the raw measurements made by the sensors to derive the raw gravity disturbance.The first step is to resample and interpret the raw measurements from the data acquisition system.Since all inertial data are recorded at a rate of 2 KHz that doesn't need for airborne gravimetry, these data are resampled by averaging over 20 samples (100 Hz).On the other hand, position and velocity of the aircraft which will be the observations of the Kalman Filter (KF), and used to compute Eötvös correction in the last step are processed with the Waypo int GrafNav software (Version 8.30) and the accelerations are obtained by double differentiation of the GPS positions.
The second step is the update of the position, velocity and attitude matrix of the aircraft from these filtered measurements.A conventional 15-state Kalman Filter is used to integrate the SINS and DGPS measurements [34].Then, the specific force in a local-level ellipsoidal frame is extracted from resampling accelerations (specific force in body frame) and the attitude matrix.
In the last step, according to the Equations ( 5) and ( 6), the raw gravity disturbance at the flight height is estimated directly from the difference between the measured specific force and the vehicle acceleration after implementing Eötvös correction and subtracting the normal gravity.
Figure 5 gives the differences between INS specific force and the derived aircraft acceleration from GPS in the vertical channel, corrected for Eötvös effect and normal gravity of the primary repeated line A and B separately.The over 10 5 mGal difference shown in the figure implies that extracting gravity disturbance from such noise data is just like looking for a needle in a haystack since the gravity disturbance does not exceed 100 mGal over distance of about 100 km [1,35].It also can be seen that the noise in some periods of the repeated line is much higher than it is in the 180° turning period of the aircraft (shown in green).This implies that the flight was implemented under the rough conditions.In fact, there were major wind direction changes especially in the region of the mouths of rivers (seen in Figure 3), which made the aircraft motion more complicated.The power spectrum of the raw gravity measurements is shown in Figure 6.As can be seen from Figure 6, much of the noise is distributed in the short wavelength of the raw data of both repeated flights.It is believed that the amplification of observation noise in the process of computing aircraft acceleration and the effects of aircraft vibration on the INS are the main part of the high frequency noise.In addition, the noise in the lower frequency, e.g., below 0.05 Hz, can be from the nature phugoid motion of the aircraft or the effect of attitude errors.It is also noted that comparing the raw data between flight A and B in both Figures 5 and 6, the noise in flight B is much stronger than flight A. This implies that extracting gravity disturbance from the raw data of flight B is more difficult than from flight A.

Test Results and Analysis
After getting the raw measurement data, EMD was employed to estimate the gravity disturbance at the flight altitude in the n-frame.Through experimentation, thirteen and fifteen IMFs were obtained for flight A and flight B, respectively.To obtain the estimation of each pure IMF, the method of using a soft thresholding was applied.The differences between two passes of each flight were used to assess internal consistency of the strapdown airborne gravity measurements.Figure 7 shows the final result of gravity disturbance denoised by EMD (shown in dash lines) and FIR (shown in thin solid lines).The blue lines denote the northbound pass and the red lines represent the southbound pass.Comparing to the raw measurement data (shown in Figure 5), the denoising effect of the used EMD is evident.However, it is well-suited to capture global trends but fails for local phenomena-some unexpected errors exist especially in flight B. Since most of the differences are distributed in the long wavelength signal and the location of the errors coincides with the mouths of several rivers around the north latitude 71.5° and 72.5° (seen in Figure 3), it is believed that the errors, which departs significantly from Gaussian white noise, are caused by aircraft motion in much dynamic environment.We have calculated the root mean square of the differences between all repeated passes by using the measure method of internal consistency proposed by Wei [7].Table 1 summarizes the estimates of the difference between two passes of each flight.The formal FIR results are also shown in Table 1 (seen in [31]).EMD Results show RMS values of 0.9 and 1.6 for flight A and B, respectively.Comparing to the FIR results (1.5 and 2.6), the corresponding improvement ratios of both profiles are more than 40%.This shows the good performance of EMD denoising.To estimate the external accuracy of the result, former airborne gravity data obtained from the LCR in 2001 and marine data around this area was used as the external reference data.Furthermore, the result of using a 160 s length FIR filter that corresponds to the same resolution as LCR data is also given.Figure 8 shows the estimated gravity disturbance for both profiles from these three types of data.Obviously, the overall trends of those three results appear to be more or less the same in both profiles.There is little difference between the EMD and FIR results.Comparing to the reference data, the result of EMD denoising shows a little smoother than FIR and is closer to the reference data.Table 2 gives the statistics of the differences as an estimation of external accuracy for those results.The difference between the RMS values for the two methods is not significant.However, for the EMD denoising, the corresponding improvement ratios of both profiles, which are about 9% and 12% respectively, at least show a good vision of using EMD to suppress the noise in airborne gravity data.

Conclusions
In this paper, we tried to apply EMD for estimating the gravity disturbance from the raw data obtained from the Greenland test of SGA-WZ.Through the computations, according to the statistics of both internal accuracy and external accuracy, although EMD denoising hasn't offered a much more accurate estimation of the gravity disturbance compared with FIR, the result at least implies EMD denoising can serve as a good tool to extract the gravity disturbance signal from noise in the strapdown airborne gravimetry.Moreover, there are still large differences at certain segments where the motion of the aircraft was much complicated.It is believed that the noise at these places departs significantly from Gaussian white noise.More sophisticated method to filter the IMFs, e.g., instead the university thresholding with a signal-adaptive one, should be focused on solving this problem in future work.

Figure 1 .
Figure 1.Main process of sifting algorithm of EMD.The procedure stops if the final residual is a monotonic function.The result is a set of IMFs and the residual.

Figure 2 .
Figure 2.An example of decomposed raw data by EMD (Only the top six IMFs are shown).The left-top plot is the first IMF, and the right-bottom plot is the sixth IMF.The noise exists in all IMFs but dominates the higher frequency IMF.

Figure 3 .
Figure 3. Flight profiles in the Greenland test and existing ground, marine and airborne coverage of the region [31].The four flight profiles shown in black (the primary flight lines are in bold) were designed in central East Greenland.Breaks in lines C and D are due to altitude changes.Other light color lines denote the marine data.The black dot represents the GPS base station and the colors show the free-air anomalies.

Figure 4 .
Figure 4.A summary of the process of using the raw measurements made by the sensors to derive the raw gravity disturbance[31].Waypoint GrafNav software is used to interpret GPS data.The 15 states in Kalman Filter are error of positions, velocities, attitudes, bias of accelerometers and gyroscopes.

Figure 5 .
Figure 5. Raw gravity measurements of repeated profile A (top) and B (bottom).The segments shown in blue and red denote the northbound and southbound separately and the green line represents the 180 degree turning time during the flights.The difference varies by over 10 5 mGal.

Figure 6 .
Figure 6.Power spectrum of the raw gravity disturbance of repeated profile A (top) and B (bottom).The noise distributes in a wide range of frequency.

Figure 7 .
Figure 7.Comparison of airborne gravity disturbance between two passes of the repeated flight A (top) and B (bottom).The blue lines denote the northbound pass, the red line denote the southbound pass.The dash lines represent EMD and the thin solid lines represent FIR.Larger differences appear around 71.5° and 72.5°.

Figure 8 .
Figure 8. Comparisons between the estimated results of FIR, EMD and reference data.The blue line denotes the FIR denoising result, the red line represents the result of EMD denoising, the green line is LCR data and the black one is marine data.(Top) is the repeated flight A and (bottom) is flight B.

Table 1 .
Statistics of the differences for internal consistency (units: mGal).

Table 2 .
Statistics of the differences for external accuracy (units: mGal).