CHAPTER 3 GLOBAL MODELS OF MEAN SEA SURFACE AND GRAVITY
3.3 Conversions from DOV to MSSH and Gravity Anomaly
We use the deflection-geoid and inverse Vening Meinesz formulae as the basic tool for computing MSSH and gravity anomaly. The deflection-geoid formula transforms DOV into geoidal height, which then yields MSSH by adding the quasi time-independent SST. Detailed derivations of these two formulae are given in Hwang (1998). These two formulae read
q
R: mean earth radius, 6371000 m is used
γ : normal gravity, based on GRS80 (Torge, 1989) H
C′and ′: Kernel functions :
, q
q η
ξ north and east components of DOV at q (dummy index)
αqp: azimuth from q to p σ : unit sphere
dσq: surface element = cosφqdφqdλ and φq,λqare latitude and longitude
The kernel functions C′andH′are functions of spherical distance only and are defined in Hwang (1998). The 1D FFT algorithm is used to rigorously implement
Equation (3.7). In the case of using a 360-degree reference field (see below), an optimal effective radius of integration in Equation (3.7) is about 110 km (about 1º at the equator). In the 1D FFT algorithm, all geoidal heights or gravity anomalies at a fixed latitude (or parallel) are computed simultaneously (Hwang, 1998), and this is why the 1D FFT algorithm is faster than the strait sum algorithm.
Because of the singularity of the kernel function C′andH′ at zero spherical distance, the inner most zone effects on geoidal height and gravity anomaly must be taken into account and are computed by
⎭⎬ s0 is the size of innermost zone, which can be estimated from the grid intervals as
π y s ∆x∆
0 = (3.9)
Formulae such as those in Equation (3.9) are based on spherical approximation.
Errors arising from spherical approximation are investigated in detail by Moritz (1980). When using the remove-restore procedure, error in using spherical approximations should be very small compared to data noises. Consider the formula of error-free LSC in the case of using ellipsoidal correction (Moritz, 1980, p. 328)
where s, l,C and sl Cll−1 are defined in Equation (3.6) and e2 is the squared eccentricity of a reference ellipsoid, which is about 0.006694 for the GRS80 ellipsoid.
The procedure in Equation (3.10) is first to remove the ellipsoidal effect of the data, l1
e2 , perform LSC computation and finally add back the ellipsoidal effect of the signal e2s1. In the case of using remove-restore procedure where a reference field is removed from the data and the residual signal is to be recovered, both l1ands1 will be very small compared to their full signals (see Moritz (1980), p. 327, where the low degree part will vanish due to the use of a reference field). For example, if the largest element (DOV) in l1 is 100 µrad, then the largest element in e2l1 will be 0.66 µrad. This value is far smaller than the noise of DOV from the multi-satellite altimetry. Furthermore, if the largest element in s1 is 100 mgal, then the largest ellipsoidal effect on gravity anomaly is 0.66 mgal, which is much smaller than the error of the recovered gravity anomaly. Of course, whether the reference field used in the remove-restore procedure will introduce additional error is another issue.
3.4 Computation and Analysis of Global MSS Model
After the tests performed in the previous sections and the selection of a set of optimal parameters, MSSH and gravity anomalies on a 2´×2´ grid were computed over the area: 80ºS-80ºN and 0º-360ºE. The computations were divided into 36 areas, each covering a 40º×40º area. A batch job was created for each of the 36 areas and this batch job creates maps of altimeter data distribution, predicted MSSH and gravity anomaly, and many statistics for detailed examinations. The final global MSSH and gravity anomaly grids are the combination of the results from the 36 areas. Figure 3.2 shows the flowchart of computation in a 40º×40º area. The most time-consuming part
in the procedure is the gridding of DOV. Only few minutes of CPU time is needed for the 1D FFT computation of geoid or gravity anomaly in a 40º×40º area on a Pentium III 600 machine. This computational procedure uses the EGM96 gravity model (Lemoine et al., 1998) to harmonic degree 360 as the reference field. In summary, reference DOV implied by EGM96 were removed from the raw DOV to yield residual DOV, which were used to compute residual geoidal heights and gravity anomalies using Equation (3.7). The final geoidal heights and gravity anomalies are obtained by adding back the EGM96-implied values.
As seen in Figure 3.2, the procedure for obtaining the global MSSH grid is more involved than that for the gravity anomaly grid. A preliminary MSSH grid was first obtained by adding the 1994 Levitus SST to the geoid grid. Because of the use of DOV as data type, it is possible that the long wavelength part of MSSH is lost in using the deflection-geoid formula. To mitigate such a loss, we first computed the differences between the along-track T/P, ERS-1 MSSH and the preliminary MSSH.
Each difference is associated with a weight, which is the inverse of noise variance.
Then smoothing and de-aliasing of the differences were made by computed the weighted median values within 15΄×15΄ cells. The weighted median values were then interpolated on a 15΄×15΄ grid using the minimum curvature method (Smith and Wessel, 1990). The final MSSH grid is obtained by summing the difference grid (which is now re-sampled into a 2΄×2΄ grid) and the preliminary grid. The resulting grids are now designated as the NCTU01 MSSH grid and the NCTU01 gravity grid.
Because of the final adjustment using T/P MSSH, The SSH from the global MSSH grid is the height above a geocentric ellipsoid with a semi-major axis=
6378136.3 m, and flattening=1/298.257222101. The geodetic coordinates of both
In addition, the normal gravity for the global gravity anomaly grid is GRS80 because we have removed the zonal spherical harmonic coefficients C20, C40, C60 and C80 of the GRS80 reference ellipsoid (cf: Torge, 1989) when computing the reference gravity anomalies from EGM96
A standard method for evaluating a MSSH grid is to compare modeled MSSH and MSSH-derived gradients with averaged SSH and gradients from repeat missions.
The first comparison is for the MSSH values from our model and from repeat missions. Table 3.2 shows the result of the comparison of MSSH values using the averaged along-track SSH from T/P and ERS-1. Also included in Table 3.2 are the comparisons for NASA/GSFC model (Wang, 2000) and the CLS model (Hernandez and Schaeffer, 2000). Compared to the NASA /GSFC and CLS models, the NCTU01 model agrees best with the T/P and the ERS-1 MSSH. The GSFC MSSH is slightly worse than the NCTU01 MSSH. On the continental shelves (depths below 200 m), all MSSH models contain large errors, which in the case of CLS has exceeded 20 cm.
Even in the median depths, the accuracy of MSSH is not very promising and is generally worse than 10 cm. All MSSH models have the best accuracy of few cm in the deep oceans. The differences in the case of ERS-1 are smaller than in the case of T/P. This is to be explained by the fact that ERS-1 has a higher data density than T/P, so the former will dominate the resulting MSSH model. As such, the MSSH model will have a better match with the ERS-1 SSH than T/P SSH. Furthermore, Figure 3.3 shows the differences between the NCTU01 and T/P MSSH. Again the differences in Figure 3.3 are large over shallow waters and small in the deep oceans. It is noted that the pattern of differences is very similar to the pattern of ocean variability in Figure 2.3.
remove DOV of
Seasat SSH Geosat SSH ERS-1 SSH ERS-2 SSH T/P SSH
north and east DOV
MSSH grid gravity anomaly grid
Figure 3.2: Flowchart for computing global MSSH and gravity anomaly grids
This shows that the MSSH is less reliably determined in areas of high ocean variability than in other areas. Such a result agrees with the expected outcome of least-squares collocation (see Equation (3.6)): data with larger noises (variability) yield less reliable results.
Table 3.2: RMS differences (in cm) between global sea surface models and T/P and ERS-1 MSSH
T/P ERS-1
0-0.2a 0.2-1 1-10 0-10 0-0.2 0.2-1 1-10 0-10
NCTU01 16.5 10.6 1.9 5.0 6.8 4.2 2.7 3.1
CLS 24.7 26.1 3.9 9.0 12.2 8.8 4.7 5.4
GSFC00 19.4 13.5 2.0 6.0 10.7 7.6 3.6 4.3
No. points 34824 23015 485762 543601 63828 48296 1517901 1630025
a Depths from 0 to 0.2 km
Figure 3.3: Difference between NCTU01 and T/P mean sea surface heights
Table 3.3 shows the RMS differences in gradients between the NCTU01, GSFC and CLS modeled gradients and the averaged along-track gradients from T/P and ERS-1 missions. The NCTU01 and GSFC models have almost the same level of accuracy, while the CLS model has a slightly poorer accuracy than the two. For all models the distribution of error is the same as that of MSSH error, namely, bad accuracy over shallow waters and good accuracy in the deep oceans. Again, for the reason of data dominance, the error in the ERS-1 case is smaller than the error in the T/P case. In the T/P case the RMS difference of 7-9 µrad has far exceeded than the standard deviation of T/P gradient, which is 1.9 µrad based on Equation (3.3) (with an averaged point spacing of 6.6 km, see Table 2.2). On the other hand, in the ERS-1 case the RMS difference of 3.2 µrad agrees well with the expected standard deviation of 4.7µrad of ERS-1 gradient.
Table 3.3: RMS differences in along-track sea surface height gradients (in µrad) derived from global sea surface models and from T/P and ERS-1 MSSH
T/P ERS-1
0-0.2a 0.2-1 1-10 0-10 0-0.2 0.2-1 1-10 0-10
NCTU01 23.3 16.7 2.4 7.1 7.1 5.0 2.8 3.2
CLS 26.4 19.9 2.9 9.0 8.7 5.3 2.8 3.3
GSFC00 22.9 16.9 2.3 7.0 7.3 4.3 2.8 3.2
No. points 33687 22814 485561 542062 61374 47071 1510808 1619253
a Depths from 0 to 0.2 km
3.5 Computation and Analysis of Global Gravity Anomaly Model
To evaluate the accuracy of the NCTU01 gravity grid, we made comparison between the predicted and shipborne gravity anomalies in 12 areas in the world oceans (see Table 3.1 for the boundaries). These shipborne gravity data are from the National Geophysical Data Center. The corresponding ship tracks are shown in Figure 3.4. The 12 areas are so selected that deep and shallow waters, and low and high latitudes are covered in the comparison. The comparing procedure and the adjustment of shipborne gravity data were detailed in Hwang et al. (1998). Briefly, long wavelength biases in the shipborne gravity anomalies were first adjusted using altimetry-derived gravity anomalies. A pointwise comparison was then made for the adjusted shipborne gravity anomalies. Table 3.4 shows the result of the comparison.
In Table 3.4, the gravity grid from Hwang et al. (1998) is also compared. It has been shown that the gravity grid of Hwang et al. (1998) has a better accuracy than that of Sandwell and Smith’s (1997). As shown in Table 3.4, the accuracy of the current gravity grid has indeed improved over the gravity grid of Hwang et al. (1998). The improvement is most dramatic in the deep oceans such as the East Pacific Rise. In the Caribbean Sea, the improvement is also significant. Other areas with significant improvements are the Ross Sea and Kerguelen Plateau, which are situated in the Antarctica. In general, large difference between predicted and shipborne gravity anomalies may be due to either or both of the following reasons: (1) bad altimeter data quality, including data scarcity, and (2) large gravity signatures, for example, in areas with trenches and seamounts. The largest difference in the Fiji Plateau and then the Mediterranean Sea should be attributed to the first reason. Due to possibly bad
tide correction and contamination of altimeter waveforms by land mass and reefs, the predicted gravity anomalies in the Fiji Plateau are not expected to be in good quality.
A recent detailed analysis of errors in altimeter-derived gravity anomalies is given by Trimmer et al. (2001). To have a dramatic improvement of gravity accuracy in such an area, a better tide model and improved determination of SSH by retracking altimeter waveforms will be needed. Improved SSHs by wave retracking have been reported in, e.g., Anzenhofer and Shum (2001), Deng et al. (2001), and Fairhead and Green (2001).
Figure 3.4: Distributions of shipborne gravity anomalies in the 12 areas where the NCTU01 gravity anomaly grid is evaluated.
Table 3.4: RMS differences (in mgal) between predicted and shipborne gravity anomalies in 12 areas
Comparison area This work Hwang et al. (1998) Improvement
(1) Alaska Abyssal 4.887 5.168 0.281
(2) East Pacific Rise 3.057 6.786 3.729
(3) Caribbean Sea 9.840 11.399 1.559
(4) Reykjanes Ridge 5.122 5.202 0.080 (5) Sierra Leone Basin 3.678 3.688 0.010 (6) Mediterranean Sea 13.026 13.480 0.454 (7) Carlsberg Ridge 7.347 7.397 0.023 (8) Kerguelen Plateau 6.017 6.931 0.914 (9) South China Sea 8.004 8.229 0.225 (10) Mariana Trench 11.561 11.814 0.253 (11) Fiji Plateau 13.365 14.148 0.783
(12) Ross Sea 7.634 8.477 0.843
CHAPTER 4
SHALLOW-WATER ALTIMETRY GRAVITY RECOVERY: DATA PROCESSING AND
COMPARISON OF METHODS
4.1 Introduction
Satellite altimetry can be useful in coastal area, such as coastal geoid modeling and offshore geophysical explorations, see, e.g., Hwang (1997), Andersen and Knudsen (2000). Recent global altimeter-derived gravity anomaly grids by, e.g., Sandwell and Smith (1997), Andersen et al. (2001), and Hwang et al. (2002), confirm that the accuracies of estimated gravity anomalies range from 3 to 14 mgals, depending on gravity roughness, data quality and areas of comparison. These papers made their comparisons with shipborne gravity anomalies mostly in the open oceans.
Research has shown that over shallow waters, altimeter data are prone to measurement errors and errors in geophysical corrections. For example, Deng et al.
(2002) shows that within about 10 km to the coastlines, altimeter waveforms are not what the onboard tracker has expected and use of ocean mode product altimetry leads to error in range measurement. The footprint of a radiometer is also large enough to make the water vapor contents measurement near shores highly inaccurate, yielding bad wet tropospheric corrections. Large tide model errors and large wave heights, among others, add to the problem of poor quality in altimeter data over shallow waters. Worst still, there is no data on land for near-shore gravity computation and this
(for example, transforming gravity anomaly to geoid requires global integration).
Thus, gravity anomaly prediction over shallow waters has been a challenging task, due mainly to problems with both data and theory. In view of these problems, this chapter compares two methods of gravity computation from altimetry discussed in (2.3.1) and (2.3.2), and a method of removing outliers in altimeter data based on a nonlinear filter is also discussed. The test area is over shallow waters in Taiwan Strait and East China sea. Figure 4.1 shows the bathymetry based on the ETOPO5 depth grid. The water west of Taiwan is a part of the east Asia continental shelf with depths less than 200 meters, while the waters east of Taiwan is deep due to the subducting Philippine Sea Plate.
4.2 Data Processing：Outlier Detection and Filtering
Outliers which are anomalous measurements in data will create a damaging result. There are methods abundant in literature for removing outliers in one-dimensional time series, see, e.g., Kaiser (1999), Gomez et al. (1999) and Pearson (2002). In this chapter we use an iterative method to remove outliers in along-track altimeter data. Consider the measured height or the differenced height as a time series with along-track distance as the independent variable. First, a filtered time series is obtained by convolving the original time series with the Gaussian function
2
where σ is the 1/6 of the given window size of convolution. The definition of Gaussian function is the same as that used in GMT (Wessel and Smith, 1995). For all data points the differences between the raw and filtered values are computed, and the
Figure 4.1: Contours of slected depths around Taiwan, unit is meter
standard deviation of the differences is found. The largest difference that also exceeds three times of the standard deviation is considered an outlier and the corresponding data value is removed from the time series. The cleaned time series is filtered again and the new differences are examined against the new standard deviation to remove a possible outlier. This process stops when no outlier is found.
We choose pass d64 of Geosat/ERM, which travels across Taiwan, to be used for the test of outliers detection. Figure 4.2 shows the ground track of pass d64 and Figure 4.3 shows the result of outlier removing. In Figure 4.3, we clearly see that erroneous differenced heights are successfully removed (discrete points) after outlier
ERS-1/gm and Geosat/gm, we do outlier detection as well as filtering for the differenced heights with a 14-km wavelength in order to reduce the data noise. There is no need of filtering for the differenced heights from the repeat missions because their data noises will be reduced due to time averaging. As an example, we choose pass a27 of Geosat/GM for testing outlier removal and filtering. Figure 4.4 shows the ground pass a27 and Figure 4.5 shows the result. Again, our algorithm removes outliers and filter the data successfully.
Furthermore, Figure 4.6 and 4.7 show the result of outlier detection along Tracks d0222 and a3083 of Geosat/GM. Both tracks pass through Peng-Hu Island in the Taiwan Strait, where most of the outliers occur and the differenced heights contain large variation. These examples from the Geosat altimeter show that differenced height is very sensitive to sudden change of height, and is particularly useful for outlier detection in the above iterative method.
We also tested this outlier detection method using altimeter data from the repeat and non-repeat missions in the East China Sea area. The altimeter data type used is also differenced height. It turns out outliers in differenced heights are very sensitive to this method, especially when along-track SSHs experience an abrupt change. A spike of SSH will translate into two large height differences. As an example, Figure 4.9 shows the result of outlier detection along two tracks of Geosat/ERM. In this case, the outliers are removed using a 28-km window size to generate the filtered time series.
Another test was carried out using data from a non-repeat mission--Geosat/GM.
Figure 4.11 shows the result of the test along two tracks of Geosat/GM. In this case, we used an 18-km window size. According to our experiments, different window sizes
for outlier detection should be used for altimeter data from different missions. Based on numerous tests, in this chapter we adopt 28 km as the optimal window sizes for the repeat missions, and 18km and 14km for non-repeat missions, respectively, over the ECS and the TS.
Figure 4.2: The ground track of pass d64 of Geosat/ERM
Figure 4.3: Raw data points (red) and outliers (blue) detected with a 28-km window for pass d64 of Geosat/ERM.
Figure 4.4: The ground track of pass a27 of Geosat/GM
Figure 4.5: Raw data points (blue) and filtered and outlier-free points (red) using a 14-km window for pass a27 of Geosat/gm.
Figure 4.6: The ground tracks of Geosat/GM d0222 and a3083
Figure 4.7: Differenced heights for passes a3083 and d0222 of Geosat/gm, crosses represent outliers.
Figure 4.8: Ascending and descending passes of Geosat/ERM for outlier detection.
Figure 4.9: Results with a 28-km window size. Outliers are not connected by the lines.
Figure 4.10: Ascending and descending passes of Geosat/GM for outlier detection.
Figure 4.11: Results with a 18-km window size. Outliers are not connected by the lines.
4.3 Results of Tests
We use multi-satellite altimeter data mentioned in chapter 2 for testing two methods (by using three kinds of altimeter data types). These data have different noise levels, which work as weights in LSC computations. In all computations, the standard remove-restore procedure is employed. In this procedure, the long wavelength part of altimeter-derived data (differenced height, height slope or DOV) implied by the EMG96 geopotential model (Lemoine et al. 1998) to degree 360 is removed. With the residual data, the residual gravity anomaly is computed and finally the gravity anomaly implied by the EMG96 model is restored. The needed isotropic covariance functions are computed using the error covariance of EGM96 and the Tscherning/Rapp Model 4 signal covariance (Tscherning and Rapp, 1974). All needed covariance functions are tabulated at an interval at of 0.01°. The local covariance values in a prediction window are scaled by the ratio between the data variance and the global variance.
For the gravity grids from methods 1 and 2, a further filtering be 2 2-D median filter improved the result. Table 4.1 and 4.2 shows the results of such filtering. Based on the testing result for method 2 in Table 4.1 and for method 1 in Table 4.2, we decide to use 16 km as the filter parameters in Taiwan Strait and to use 20 km as the filter parameters in East China Sea.
To evaluate the accuracy of the gravity anomalies derived from three different
To evaluate the accuracy of the gravity anomalies derived from three different