Experimental Investigations on Airborne Gravimetry Based on Compressed Sensing

2.1. Airborne Gravimetry Technology

An object's gravity is the composition of forces of gravitation caused by the Earth and other celestial bodies and the inertial centrifugal force caused by the Earth's rotation. The non-uniform distribution of the density of the Earth's interior matter makes the gravity vary with the position. In gravity prospecting, the gravity variations caused by the non-uniform density distribution of the Earth's interior rocks and minerals are called gravity anomalies. In fact, the airborne gravity measurements include two parts: the gravity anomaly (denoted as the free air anomaly here) and the normal gravity (which is the reference gravity field of a conventional ellipsoid). Therefore, the gravity anomaly can be expressed as:

The principle of strap-down airborne scalar gravimetry is based on Newton's equation of motion in the gravitational field of the Earth, utilizing the principle of relative gravity measurement. First, we must use a terrestrial gravimeter to connect national gravity and a point on the parking apron to obtain its gravity and take it as the gravity reference point. Before the aircraft takes off, static gravimeter data must be recorded so that gravity observations in the air can be associated with the gravity reference point on the parking apron. Then, the strap-down airborne scalar gravimetry model can be written as:

Obviously, Equation (2) is an equivalent form of Equation (1), where Δg is the gravity anomaly to be determined, ν̇_D is the vertical component of the vehicle acceleration obtained from GPS, f_D is the vertical component of the specific force measured by the accelerometers of an inertial measuring unit, $f_D^0$ is the vertical component of the specific force obtained from the static data on the parking apron, called base reading, g_b is the gravity of reference point on the parking apron and Δa_E includes all kinds of the error correction [16].

2.2. Mathematical Model of CS

Compressed sensing (CS) is an exciting, rapidly growing field that has attracted considerable attention in fields as diverse as electrical engineering, applied mathematics, statistics, sensor technology and computer science. CS offers a framework for simultaneous sensing and compression of finite dimensional signals. Quite surprisingly, it predicts that sparse high-dimensional signals can be recovered from highly incomplete measurements by using efficient algorithms. CS also holds promise for increasing resolution by exploiting the signal structure. Especially, reducing the sampling rate or increasing resolution in airborne gravimetry can improve survey efficiency, increase data transfer rate and improve data quality.

Let z be an unknown original gravity anomaly in Rⁿ (which can represent a 1-D or 2-D gravity anomaly of interest). Suppose that we have m linear measurements of z with the form:

When Ψ is invertible, CS exploits the sampling gravity anomaly data reconstruction by solving a problem of the form:

2.3. OMP Algorithm

Let A = ΦΨ⁻¹ ∈ R^m×n; Equation (5) can be rewritten as:

Equation (6) is the zero-norm sparse optimization problem and can be solved with greedy algorithms. Greedy algorithms rely on iterative approximation of the signal coefficients and support, either by iteratively identifying the support of the signal until a convergence criterion is met, or alternatively, by obtaining an improved estimate of the sparse signal at each iteration that attempts to account for the mismatch to the measured data. Greedy algorithms are classified as greedy pursuits and thresholding type algorithms. Greedy pursuit algorithms include Matching Pursuit [18], OMP [19], Conjugate Gradient Pursuit [20], Stagewise Orthogonal Matching Pursuit [21] and Regularized OMP [22,23], etc.

OMP is a greedy pursuit algorithm. Given that z is k-sparse in the transform domain Ψ, i.e., x has at most k nonzero elements and A = ΦΨ⁻¹ = [a₁,a₂,⋯,a_n] whose columns denoted by a_i ∈ R^m, the vector y = Ax is a linear combination of k columns from the matrix A. In the language of sparse approximation, we say that y has an k-term representation over the dictionary A. Therefore, OMP can be used for recovering sparse signals.

To identify the signal x, we need to determine which columns of A participate in the vector y. The idea behind the OMP algorithm is to pick columns in a greedy fashion. At each iteration, we choose the column of A that is most strongly correlated with the remaining part of y. Then we subtract off its contribution to y and iterate on the residual. One hopes that, after k iterations, the OMP algorithm will have identified the correct set of columns. The major advantages of OMP are its speed and its ease of implementation. The procedure of OMP can be carried out by means of the basic steps as follows [15]:

Given the matrix A, the vector y and the sparsity level k of the signal x:

• Step 1: Initialize the residual r₀ = y, the index set Λ₀ = Ø and the iteration counter t = 1.

• Step 2: Find the index λ_t that solves the easy optimization problem:

  $${\lambda }_t=arg{max}_{i=1,\cdots ,n}\left|\langle r_{t-1},a_i\rangle \right|$$

  (7)

  If the maximum occurs for multiple indices, break the tie deterministically.

• Step 3: Augment the index set and the matrix of chosen atoms: Λ_t = Λ_t−1∪{λ_t} and Ω_t = [Ω_t−1 a_λt]. We use the convention that Ω₀ is an empty matrix.

• Step 4: Solve a least squares problem to obtain a new signal estimate:

  $$s_t=arg{min}_s{\Vert y-{\Omega }_{t}s\Vert }_2={({\Omega }_t^T{\Omega }_t)}^{-1}{\Omega }_t^{T}y$$

  (8)

• Step 5: Calculate the new residual:

  $$r_t=y-{\Omega }_t{s}_t$$

  (9)

  It is important to note that the residual r_t is always orthogonal to the columns of Ω_t. Provided that the residual r_t−1 is nonzero, the algorithm selects a new atom at iteration t and the matrix Ω_t has full column rank. For such a case, the solution s_t to the least squares problem in Step 4 is unique.

• Step 6: Increment t, and return to Step 2 if t < m.

• Step 7: The estimate x̂ for the signal x has nonzero indices at the components listed in Λ_m. The value of the estimate x̂ in component λ_t equals the tth component of s_m.

• Step 8: The estimate for the unknown signal z is: ẑ = Ψ⁻¹x̂.

3.1. Flight Test and Data Preprocessing

The strap-down airborne scalar gravimeter called SGA-WZ mentioned in this paper is the first system with this type in China. It was developed by the Laboratory of Inertial Technology of the National University of Defense Technology [13]. This system consists of a high-performance strap-down inertial navigation system (SINS), a GNSS receiver, an anti-vibration system, a data logger and post-processing software. The major advantage of this system is its reliability and robustness in operations. A photograph of the system is shown in Figure 1.

The flight tests were carried out in Shandong Province using SGA-WZ from April 2010 to May 2010. The hardware was installed in the aircraft six days before the flight tests. Two GNSS ground stations were located near the airport where the aircraft took off and landed. GNSS receivers installed on the ground and aircraft can be used to determine the vehicle position, velocity, and acceleration. The strap-down airborne scalar gravimeter onboard a Cessna 208 aircraft was used to collect the data. Figure 2 shows the Cessna 208 aircraft, which is a fixed-wing small aircraft.

The pilots controlled the aircraft with the autopilot, and the test was implemented in days with good weather to minimize the effects of air turbulence. The average flight altitude was approximately 400 m above sea level with a fluctuation of 20 m. The average speed during the flight was 60 m/s. The sampling rate of raw SINS readings was 100 Hz and 2 Hz for the GNSS sampling rate. After being installed in the aircraft, SGA-WZ worked all day for over a month. In the whole flight test campaign there were eight flights. The first and second flights were repeated lines that flew along the same trajectory to test the repeatability of the system, and the valid length of each repeated profile was about 100 km. The other six flights left were grid flights consisting of three flights of survey lines and three flights of control lines. The spacing between survey lines was about 2 km and the spacing between control lines was about 9 km. Figure 3 shows the grid flight lines used in the test [16].

As previously mentioned, the gravity anomaly was determined from the difference between the specific force and the vehicle acceleration using Equation (2). Because of the phugoid modes of the aircraft and the atmospheric turbulence, there were many high frequency noises in the gravity anomaly. The noises should be eliminated by the FIR low-pass filter. When the velocity of aircraft is fixed, the longer the filter length, the more high frequency noises can be eliminated, but the spatial resolution will be lowered. The test results have revealed that the internal accord precision is 1.50 mGal with a spatial resolution of 6 km. Figure 4 shows the gravity anomaly of line 1 in the first flight before low-pass filtering. Figure 5 shows the gravity anomaly after 160 s cutoff low-pass filtering.

3.2. Sparsity Analysis

The first and second flights were repeated lines surveys denoted by F401 and F402, respectively. Each of them consisted of six lines. Figure 6 shows the gravity anomaly curves of F401 and Figure 7 shows the gravity anomaly curves of F402. In the two flights, the repeated lines were surveyed along an east-west direction.

Based on the theory of CS, this study analyzes the sparsity of airborne gravimetry by DFT. Figure 8 shows the amplitude spectrum of all lines in F401 and Figure 9 shows the amplitude spectrum of all lines in F402. From the two figures, we can see that the gravity anomaly is mostly distributed in the low-frequency domain and it has few high-frequency components. The results thus reveal that the gravity anomaly is sparse in the DFT domain, which satisfies the precondition for application of CS.

3.3. OMP Reconstruction Results

To evaluate the OMP algorithm performance and the gravity anomaly data reconstruction precision, we now define the performance indices with standard deviation and SNR. SNR can be calculated by:

In this study, the compressed measurement number was set as m = 432, i.e., the compressed sampling rate was approximately 14%. Figure 10 shows the reconstruction result of line 1 in F401 and the reconstruction error is shown in Figure 11. In Figure 10, the red dashed line denotes the reconstruction gravity anomaly and the blue solid one denotes the original gravity anomaly. The two figures show that the OMP results approximate the original gravity anomaly very well and the reconstruction error is very small. However, we also note that a boundary effect exists in the result and this problem should be solved in the future.

For comparison purposes, we also made the gravity anomaly data reconstruction with NIPM, which is one of the commonly used methodologies for data reconstruction. Figure 12 shows the reconstruction result of line 1 in F401 and the reconstruction error is shown in Figure 13. The comparison of the algorithm performance between OMP and NIPM is shown in Table 1. It shows that the each performance index from OMP is superior to the corresponding one from NIPM. The standard deviation of the reconstruction error and SNR of OMP are 0.03 mGal and 56.48 dB, respectively. Correspondingly, the standard deviation of reconstruction error and SNR of NIPM are 0.15 mGal and 42.29 dB, respectively.