Chapter 2 Data processing in scalar airborne gravimetry: theory and technique
2.6 Summary and discussion
12
yy xx
xy
xy S S
K S
, with 0Kxy 1 (2.46) where is frequency, Sxy the cross spectral density between x and y, and Sxx
and S the auto spectral density of x and y, respectively. yy
The spectral coherence is a statistic that can be used to estimate the power transfer between x and y. The coherence values will satisfy 0Kxy 1. If x and y are completely unrelated the coherence will be zero. On the contrast, it will be one. If it is less than one but greater than zero it indicate either noise is entering the measurements or the relation between x and y is not linear. In the case that Kxy 0.5, x and y are said to be coherent at a particular frequency.
2.6 Summary and discussion
This chapter presents the basic principles, data processing and methods of accuracy analysis of scalar airborne gravimetry. A theory, with formulae based on the global average power of gravity anomaly, is presented and can be used to predict a theoretical resolvable wavelength at a given flight altitude. Focusing on the GPS-related error sources, the required accuracy for a one-mgal accuracy of airborne
37
gravity value was derived for a case in Taiwan. All the programs developed by the space geodesy laboratory, department of civil engineering, NCTU for the airborne gravity data reduction have been tested and produced reliable results using data from airborne gravity surveys. The flight designs and results of the airborne gravity surveys can be found in Chapter 3. One module that may be improved is filter technique. An
“optimum” filter will be one that can minimize data noises while preserving the finest spatial resolution. Our program can easily adopt a new filter if necessary.
38
39
Chapter 3
Multiple-altitude airborne gravity surveys
Taiwan’s terrain is complex and mostly inaccessible for gravity survey. Over 75% of Taiwan’s terrain is covered with hills and high mountains, with the highest point being nearly 4000 m. The variations of gravity anomaly around Taiwan are quite large, varying from -250 mgal over the trench east of Taiwan to 450 mgal over the Central Range of Taiwan. The old surface gravity data are sparsely distributed and uncertainties in the gravity datum and the coordinate system (Hwang, 1997) (see Figure 3.1). For such applications as geoid modeling, vertical datum determination, geological study and ocean current determination, a dense, accurate gravity data set is indispensable. To this end, MOI, Taiwan, sponsored three airborne gravity campaigns over: (1) the Taiwan Island, (2) the Kuroshio Current and (3) the Taiwan Strait and Dongsha Atoll, as parts of NFS from 2004 to 2009. The field work and data processing were mainly carried out by NCTU, with beginner’s training obtained from DNSC. The survey area covers the Taiwan Island, the surrounding waters and an island in the South China Sea.
This chapter gives an overview of these airborne gravity surveys in Taiwan. The methodology for the data processing of airborne gravimetry in Chapter 2 is applied to the three airborne gravity surveys. These airborne gravity data are then used for the computation of Taiwan geoid model in Chapter 5.
40
Figure 3.1: Distributions of the terrestrial and shipborne gravity data before 2004.
3.1 Description of the airborne gravity surveys
The airborne gravity surveys in Taiwan were all made using a L&R Air-Sea Gravity System II gravimeter (serial number: S-133) (LCR, 2003). Table (3.1) summarizes the specifications of the L&R gravimeter. The King-Air Beechcraft-200 and Beechcraft-350 aircrafts were modified to accommodate the gravimeter and a Trimble 5700 GPS receiver with gravity and GPS data sampled at 1-Hz. The aircraft trajectories were determined by means of the GPS kinematic positioning technique with several ground GPS tracking stations and the GPS software, Bernese 5.0 (Beutler
41
et al., 2004; Rolf et. al., 2007). Both the gravimeter and the aircrafts are owned by MOI, Taiwan.
Furthermore, the decision of the survey plan was based on the allocated time of the aircrafts (support of other missions), Taiwan’s topography (highest point 3952 m), project duration, and a political factor (not to over the center of the Taiwan Strait).
Table 3.1: Specifications of the L&R Air-Sea Gravity System II.
Resolution 0.01 mGal
Accuracy <1.00 mGal
Size 71 x 56 x 84 cm
Weight Meter: 86 kg, UPS: 36 kg
Power requirements
240 watts(avg), 450 watts(max), 80-265 VAC, 47-63 Hz
Recording rate 1-Hz
3.1.1 Taiwan Island airborne gravity survey at the altitude of 5000 m
The Taiwan Island airborne gravity survey was carried out at an altitude of 5000 m over the period of May 2004 to March 2005. The goal of this project is to obtain evenly distributed free-air gravity data covering entire Taiwan Island, especially over the Central Range of Taiwan.
The field work lasted 43 days, including 3 days of re-flights, and exceeded 200 hours. Since Taiwan is long in the north-south direction and short in the west-east direction, most of the flight lines are in the north-south direction to avoid excessive turns of flight and to reduce unusable data in the beginning of flight line. Besides,
42
north-south lines contain the minimal Eötvös effects. In general, north-south lines yield better gravity accuracy than that from west-east lines.
Figure 3.2: Numbers of survey lines and GPS tracking stations in Taiwan Island airborne gravity survey (★: airport; ▲: GPS tracking station).
Based on this consideration, the numbers of survey lines in the north-south, east-west, northeast-southwest and northwest-southeast directions are 64, 22, 10 and 6 respectively, see Figure (3.2). The cross-line spacing is 4.5 km for all survey lines, except the east-west lines, which are spaced at 20 km. The west-east lines are mainly
43
for crossover analysis. The speed of flight is 306 km/hour relative to the surface. The gravity readings sampled at a 1-Hz correspond to an 85-m sampling interval on the surface. For the kinematic positioning of the aircraft, GPS data at 7 permanent GPS tracking stations around Taiwan and 1 tracking station (SNAM, Figure 3.2) near the Taichung airport were collected. The gravity value at the aircraft parking spot was determined using relative gravity measurements by a Graviton-EG gravimeter (LCR, 2002) and tied to an absolute gravity value at the FG5 gravity station of Taichung. The standard error of this absolute gravity value is 0.04 mgal based on a relative gravity network adjustment.
3.1.2 Kuroshio Current airborne gravity survey at the altitude of 1500 m
The Kuroshio Current airborne gravity survey was flown over the period of March 2006 to August 2008. This project aims to map the Kuroshio Current east of Taiwan by the combination of airborne gravity data, surface gravity data and satellite altimetry data. The gravity data are used to compute an accurate geoid model over the Kuroshio area east of Taiwan. Satellite altimetry data to be used are from the Geosat, ERS-2, ENVISAT, T/P and JASON-1 missions. The altimetry data are improved using waveform retracking, filtering and outlier rejection techniques. With a high resolution and accurate geoid model, a two-dimensional (2D) mean velocity field of the Kuroshio east of Taiwan may be computed. The east and north components of current velocity at a crossover point of two satellite ground tracks can be computed for all altimetry missions. Time series of current velocity at crossover points are used to see temporal variations of currents. The crossover points offer opportunity to investigate phenomena such as the intrusion of the Kuroshio to the South China Sea, the intrusion of the Kuroshio to the East China Sea and the branching out of the Kuroshio to the
44
east before reaching the I-Lan Ridge. The detailed description can be found in Chapter 6.
Figure 3.3: Numbers of survey lines and GPS tracking stations of Kuroshio Current airborne gravity survey (★: airport; ▲: GPS tracking station).
The Kuroshio Current survey contains 36 and 7 lines in the north-south and east-west directions, respectively (see Figure 3.3). The cross-line spacing is 5 km for all north-south lines and 60 km for east-west lines. All the lines are distributed over the Kuroshio Current east of Taiwan. To increase the resolution of gravity data, the
45
flight altitude is set to a 1500 m. The speed of flight is 280 km/hour with 1-Hz sampling rate, corresponding to a 77 m interval on the surface. A kinematic GPS network with 1 roving station and 6 fixed GPS tracking stations is used to determine the position of the aircraft. The absolute gravity value at airport is measured in the same way as in the Taiwan Island airborne gravity survey.
3.1.3 Taiwan Strait and Dongsha Atoll airborne gravity survey at the altitude of 1500 m
The Taiwan Strait and Dongsha Atoll is a typical shallow-water area, defined as waters with a depth < 500 m. In this area the shipborne gravity data were sparsely distributed, and were collected by the research vessels studying marine geophysics around Taiwan (Hwang and Hsu, 2008). Altimeter-derived gravity is the most important source in this area. However, altimeter-derived gravity in shallow waters can be seriously degraded due to bad tidal correction, bad wet tropospheric correction, large sea surface variability and contaminated altimeter waveforms (Hwang and Hsu, 2004; Deng, 2003). Therefore, it is expected the gravity data here can be improved by airborne gravimetry. This motivates the Taiwan Strait and Dongsha Atoll airborne gravity survey.
In 115 flight hours over 2008-2009, the gravity values at 1500 m above sea level over the Taiwan Strait and Dongsha Atoll were collected using the same gravimeter aircrafts in the other two surveys. The spacing along the 54 north-south and 15 west-east lines are 5 and 25 km, respectively. The flight speed and data interval are 280 km/hour and 77 m, which are the same as those for the Kuroshio airborne survey.
The 7 GPS ground tracking stations scattered at western coasts of Taiwan were used to determine the aircraft positions. Figure (3.4) and (3.5) show the flight plan of the Taiwan Strait and Dongsha Atoll airborne gravity survey.
46
Figure 3.4: Numbers of survey lines and GPS tracking stations of West Taiwan airborne gravity survey (★: airport; ▲: GPS tracking station).
47
Figure 3.5: Numbers of survey lines of Dongsha Atoll airborne gravity survey and the geographical location associated with Taiwan Island (▲: GPS tracking station).
3.2 Accuracy assessment of position, velocity and acceleration
GPS not only provides the precise position of aircraft but also the velocity and acceleration for data processing. Here, accuracy assessments of position, velocity and acceleration using the GPS data collected in the multiple altitudes airborne gravity surveys are presented. The positions of aircraft were determined by the Bernese 5.0 software (Beutler et al., 2004; Rolf et. al., 2007). Ten sessions of GPS data from the Kuroshio and Taiwan Strait surveys at 1500 m, and another ten from the Taiwan Island survey at 5000 m were selected for the computation of position, velocity and acceleration and for accuracy assessments. Each session of about 4 hours in data length was divided into two independently processed sub-sessions with a 30-minute overlap. The 20 overlapping sessions were selected so that the tracks of the overlaps cover coastal plains, high mountains, and oceans.
Figure (3.6) shows a typical example of position differences in an overlapping
48
session, which vary slowly and contain spikes and discontinuities. Figure (3.7) shows the PSD of position differences for the three components. The spectral index, i.e., the value of the exponentxinS(f) f x , where S is PSD and f is frequency, is -0.7,-0.8 and -1.4 for the north, east, and vertical components, respectively. Therefore, the PSD approximately indicates flicker noise (x1) for the north and east components, and a combination of flicker noise and a random-walk process (x2) for the vertical component. Thus, the overlapping differences of position estimates are dominated by the long-wavelength errors in GPS positioning. These results are consistent with other studies, e.g., Verdun et al. (2003) represent the combined effect of long and short wavelength GPS errors. Furthermore, values of the PSD of the vertical component are 0.5 to 1.5 orders of magnitude smaller than the PSD of the north component, but similar to the PSD of the east component at high frequencies.
Values of the PSD of the north component are 2 to 4 smaller than the east component at all frequencies.
The overlapping difference can be regarded as the internal accuracy of GPS positioning. It was found that spikes and discontinuities were associated with phase ambiguity changes, i.e., change of visible satellites and/or cycle slips in one or more of the baselines. The position differences in Figure (3.6) are mainly caused by the differences in the estimated common parameters associated with the overlapping sub-sessions, e.g., phase ambiguities and tropospheric parameters. Since these common parameters will remain unchanged within a certain time span, their effects will be reduced upon differentiation. That is, we would expect that the differences in position-derived velocity and acceleration (by numerical differentiation) are less affected by the differences in common parameters.
49
Figure 3.6: Differences of position at an overlapping session in the north, east, and vertical components. The left vertical scale is for the north and east components, and the right vertical scale for the vertical component. Discontinuities occur at seconds 615 and 699 (Hwang et al., 2007).
Figure 3.7: Power spectral densities of the differences (Hwang et al., 2007).
50
Table 3.2: Averaged RMS differences of 10 overlapping sessions of GPS kinematic solution from the Kuroshio and Taiwan Strait surveys at the altitude of 1500 m.
North East Vertical
Position (m) 0.137 0.268 0.421
Velocity (m·s-1) 0.0064 0.0081 0.0108
Acceleration (mgal) 93.88 74.53 906.28
Table 3.3: Averaged RMS differences of 10 overlapping sessions of GPS kinematic solution from the Taiwan Island survey at the altitude of 5000 m.
North East Vertical
Position (m) 0.073 0.208 0.293
Velocity (m·s-1) 0.0003 0.0003 0.0013
Acceleration (mgal) 23.56 26.56 104.85
Table (3.2) and (3.3) lists the average RMS differences in position, velocity, and acceleration from the 10 overlapping sessions at 1500 m and 5000 m, respectively.
The overall positioning accuracy is of the order of decimeter, with the vertical component being the largest. The RMS velocities show that the long wavelength positioning errors does not propagate to the errors in velocity and acceleration. The statistics of differences in velocity and acceleration at the altitude of 1500 m are larger than those at the altitude of 5000 m. The reasons are (1) the measuring environment at lower altitude contains more turbulence and (2) the GPS signals at 1500 m contain more atmospheric influence such as tropospheric errors than the case at 5000 m.
Assume that the RMS differences in Table (3.3) are equal to the standard errors of the respective quantities. In theory, a standard error of 0.293 m in vertical position
51
0 m·s-2 in vertical acceleration. However, the standard errors in Table (3.3) are much smaller than these two numbers. This implies that differentiation of positions does reduce the effects of long-wavelength errors of positions on velocity and acceleration (Verdun et al., 2003).
Table 3.4: Standard deviations of differences at an overlapping session of 5000 m in the north and east velocity components and the vertical acceleration with different filter window size (Hwang et al., 2007).
Window
120 5000 0.000030 0.000042 0.000015
180 7500 0.000027 0.000037 0.000012
240 10000 0.000026 0.000034 0.000010
Furthermore, to see the performance of the Gaussian filter and the effect of filter window size on the GPS-derived velocities and accelerations we performed the following tests using the 20 overlapping sessions. The results in Table (3.4) are based on an overlapping session. We assume that the aircraft velocity is 300 km/hour, and the equivalent spatial resolution of airborne gravity data is just the multiplication of the aircraft velocity and the half-window size. In general, with a window size greater than 60 seconds, the standard deviations of velocities are at the 0.005 m·s-1 level and
52
the standard deviations of vertical acceleration are at the mgal level. A window size larger than 200 seconds yields a sub-mgal standard deviation. Through Table (3.4) we can initially understand how to select an adequate window size for the filtering of airborne data. In an actual flight, the aircraft may experience turbulence and GPS data may contain cycle slips, so the window size to be used will vary from one case to another.
3.3 Detection and down-weighting of outliers
Czombo (1994) discussed the impact related to the GPS errors on the filtered airborne gravity. The position accuracy of 0.5 m will be satisfactory for the reduction of airborne data. The reality is that it’s not possible to expect a smooth flight all the time but slight turbulent. The cycle slips from GPS and sudden position jumps caused by the atmospheric interference may produce erroneous data within an entire survey line. Inferior or erroneous data will lead to the results containing artifacts and in turn false interpretations of the underlying geophysical phenomena. It is necessary to detect and remove these errors prior to the filtering.
In statistics, an outlier is that an observation is numerically distant from the rest of the data (Barnett and Lewis, 1994). Outliers can occur in any distribution, but they often indicate measurement errors or that the data distribution has high spike. In the former case the statistics that are robust to outliers are often used to reduce the influence of outliers, while in the later case a hypothesis of normal distribution for the data is used. The 3-sigma rule is a simple approach to determine the gross errors from the observations
The raw data, containing high frequency noise, were first filtered using a twice filter window wider than that used to produce the final results. The local differences (within a filter window width) between raw and filtered data were used to compute a
53
standard deviation. By the 3-sigmal rule, each observation is examined to be an outlier or not. Figure (3.8a) shows the raw gravity anomalies along one profile, Line 17, form the Kuroshio airborne gravity survey. The outliers along Line 17 are shown in Figure (3.8b).
However, a rough down-weighting for an observation will not only affect a number of observations before and after the outliers but also probably destroy the time series of data. Alberts (2009) suggested applying a robust estimation approach to the outliers, such as M-estimation (Huber, 1981). In this study, the new weight for an outlier is reduced using the “Danish Method” as (Leick, 2004)
local differences and the difference between the raw and filtered value. An outlier with a larger difference gets a smaller weight.
After that, a Gaussian filter is used for the new weights to obtain the filtered weights. A filter window about 10-20 second is suitable, depending the duration of turbulence in the flight. Figure (3.9) shows the weights of all observations along Line 17 after down-weighting processing and after filtering. A corrected gravity data can be obtained by multiplying the filtered weights and raw gravity data. The differences between raw and corrected gravity data and the corrected gravity data are shown in Figure (3.10).
54
Figure 3.8: (a) Raw gravity anomalies along Line 17 from the Kuroshio airborne gravity survey (b) the outliers identified by the 3-sigmal rule.
55
Figure 3.9: (a) New weights of all observations after the down-weighting processing;
(b) filtered weights of all observations using a filter width of 15 s.
56
Figure 3.10: (a) The differences between raw and corrected gravity data; (b) the corrected gravity data.
57
Table 3.5: Statistics of differences between the upward-continued surface gravity and the airborne gravity filtered by the Gaussian, Butterworth and HHT filters along Line 17 from the Kuroshio survey at 1500 m (unit: mgal).
Filter type Min Max Mean STD RMS
Gaussian -5.99 11.84 0.10 3.56 3.56
Butterworth -5.72 10.33 0.40 3.67 3.70
HHT -17.91 15.14 0.08 4.96 4.96
Gaussian with outlier detection
-7.31 10.92 0.03 2.87 2.87
3.4 Filter evaluation
Two low-pass filters, Gaussian and Butterworth, and the new method, HHT, were tested. The upward-continued surface gravity data are computed for the comparison of the filtered results. It is decided to adopt a filter width of 150 s for the low-pass filter. This corresponds to a 6-km spatial resolution (half-wavelength) at the altitude of 1500 m. Figure (3.11) shows the 3nd-iterative filtered gravity using the Gaussian, Butterworth and HHT filters from the corrected gravity data (Figure 3.10b).
The statistics of differences between the upward-continued gravity and the filtered ones are summarized in Table (3.5). Before the outlier detection, the best result of 3.56 mgal was obtained with the Gaussian filter.
In Figure (3.11) the turbulences in the flight make undesired results around 750 s and 2200 s. This is due to turbulence that generates large vertical accelerations that
In Figure (3.11) the turbulences in the flight make undesired results around 750 s and 2200 s. This is due to turbulence that generates large vertical accelerations that