Improvement of the geoid model over Japan using integral formulae and combination of GGMs

Abstract

An improved high-resolution gravimetric geoid model covering the four main islands of Japan (Hokkaido, Honshu, Shikoku, and Kyushu) was developed on a 1 × 1.5 arc-minute grid from EGM2008, GOCO02S/EGM2008, and terrestrial gravity data. A modified form of the Stokes-Helmert scheme was applied for the determination of the geoid using an empirically determined optimal spherical cap size. Handling of the topographical effects on gravity was accomplished using the integral formulae of Martinec and Vaníček, which were found to be more suitable for geoid modeling over Japan than the classical formulae. EGM2008 was used in Hokkaido, Honshu, and Kyushu Islands, whereas a combination of GOCO02S and EGM2008 was used in Shikoku Island and its immediate surroundings. The global geopotential models (GGMs) used in this study were chosen based on our earlier evaluation of the performance of EGM2008- and GOCE-related GGMs in Japan. In comparison with the previous geoid model for Japan, our new model shows an improvement in the standard deviation from ±8.3 to ±7.5 cm.

Background

Geoid determination over Japan has been the focus of numerous studies over the past few decades (e.g., Ganeko 1976; Kuroishi 1995, 2001a, b,2009; Fukuda et al. 1997; Kuroda et al. 1997; Kuroishi et al. 2002; Kuroishi and Denker 2001; Kuroishi and Keller 2005; Odera et al. 2012; Odera and Fukuda 2013), although precise geoid modeling remains a challenge. In this study, an improved high-resolution gravimetric geoid model is developed from EGM2008 (Pavlis et al. 2012), GOCO02S (Goiginger et al. 2011) combined with EGM2008 (GOCO02S/EGM2008), and local terrestrial gravity data. EGM2008 is complete to spherical harmonic degree and order 2,159 and contains additional coefficients extending to degree 2,190 and order 2,159.

The choice of global geopotential models (GGMs) used in this work is based on our earlier evaluation of GGMs that have recently been developed over Japan (Odera and Fukuda 2013). In that study, the performance of EGM2008 and the Gravity field and steady-state Ocean Circulation Explorer (GOCE)-related GGMs that were recently released by the European Space Agency were assessed, the latter being based on 2, 8, and 12 months of data collected by GOCE. This assessment was based on geoid undulations, obtained from 816 GPS/leveling points, and free-air gravity anomalies, obtained from 6,951 first-order gravity points. The distribution of these data over the four main islands of Japan is shown in Figure 1.

Figure 1

GPS/leveling (big black dots) and first-order gravity (small red dots) over Japan.

The evaluated GOCE and related GGMs included GOCE-DIR1, 2, 3 (Pail et al. 2011; Bruinsma et al. 2010), GOCE-TIM1, 2, 3 (Pail et al. 2011; Pail et al. 2010b), GOCE-SPW1, 2 (Migliaccio et al. 2011), and GOCO01S, 02S (Goiginger et al. 2011; Pail et al. 2010a). The performance of EGM2008 and GOCE-related GGMs over Japan was shown by Odera and Fukuda (2013) to be almost identical at a harmonic expansion to degree 150, although EGM2008 performs better at the 180-, 210-, and 240-degree spectral bands. Furthermore, comparisons over the four main islands reveal that EGM2008 performs better than the GOCE and related GGMs in Hokkaido, Honshu, and Kyushu at 180, 210, and 240 degrees. However, the GOCE and related GGMs perform better than EGM2008 in Shikoku.

Odera and Fukuda (2013) developed a gravimetric geoid model for Shikoku and the surrounding area from GOCE-TIM3/EGM2008 and terrestrial gravity data, with the aim of determining the magnitude of improvement in the model by using GOCE-only data. However, GOCO02S/EGM2008 performs slightly better than GOCE-TIM3/EGM2008 in this region. Therefore, in this work, a combination of GOCO02S (up to 180 degrees) and EGM2008 (from 181 to 2,190 degrees) is used to improve the geoid model over the Shikoku area, and EGM2008 is used in Hokkaido, Honshu, and Kyushu. Ship-track gravity data and satellite altimetry-derived marine gravity anomalies in the coastal areas were not considered because of the use of a high-resolution gravitational model (EGM2008). Furthermore, the interest of this study is limited to the land areas.

The classical Moritz formula (Moritz 1980) and a planar formula (Wichiencharoen 1982b) were used to compute the direct topographical effect (DTE) and the primary indirect topographical effect (PITE), respectively, in the previous geoid model for Japan. These two formulae (classical Moritz and planar) are hereafter referred to as the classical formulae and are given respectively as

$$\mathit{C}=\frac{\mathit{G}{\mathit{\rho }}_0{\mathit{R}}^2}{2}{\displaystyle \underset{\mathit{\sigma }}{\iint }\frac{{\left(\mathit{H}-{\mathit{H}}_{\mathit{P}}\right)}^2}{{\mathit{l}}_{\mathit{o}}^3}}\mathit{d}\mathit{\sigma },$$

(1)

and

$${\mathit{N}}_{\mathrm{\text{ind}}}=\frac{-\mathit{\pi }\mathit{G}{\mathit{\rho }}_0{\mathit{H}}_{\mathit{P}}^2}{\mathit{\gamma }}-\frac{\mathit{G}{\mathit{\rho }}_0{\mathit{R}}^2}{6\mathit{\gamma }}{\displaystyle \underset{\mathit{\sigma }}{\iint }\frac{{\mathit{H}}^3-{\mathit{H}}_{\mathit{P}}^3}{{\mathit{l}}_{\mathit{o}}^3}}\mathit{d}\mathit{\sigma },$$

(2)

where C is the DTE, N _ind is the PITE, G is the Newtonian gravitational constant, ρ₀ is the constant topographic density, H _P is the orthometric height of the computation point, H is the height of the running point, σ is the surface integration element, and l _o is the horizontal distance between the computation point and the running point.

The current geoid model for Japan has been developed using the integral formulae of Martinec and Vaníček (1994a, b) for the computation of DTE and PITE. Details concerning the procedure and data sets used to develop the previous geoid model can be found in Odera et al. (2012), these being generally the same as for the current computations, with the exception of topographical reductions.

This study further assesses the appropriate topographical gravity reduction procedure for precise geoid modeling in Japan. Gridding of the residual gravity anomalies is accomplished by kriging (Krige 1951) on a 1 × 1.5 arc-minute grid, Stokes's integral formula (Stokes 1849; Heiskanen and Moritz 1967) is used for geoid determination, and the modified Stokes's kernel proposed by Meissl (1971) is applied to minimize truncation errors.

Methods

Gravity data reductions and gridding

Most of the gravity data used in this study was obtained from a database developed by Nagoya University and other organizations that covers the southwestern parts of Japan (e.g., Shichi and Yamamoto 2001a, b). Another set of gravity data, mostly along the main leveling network covering the four main islands (Hokkaido, Honshu, Kyushu, and Shikoku), was provided by the Geospatial Information Authority of Japan (GSI). In total, there are 98,670 observed gravity data points in the study area that were obtained during the observation period approximately 1910 to 2001 (Figure 2). It should be noted that the northern parts of Japan (above latitude 37°N) are sparsely populated with gravity data points, whereas the southern parts are densely populated; however, most of these data was clustered together and even the well-covered areas contain wide data gaps.

Figure 2

Distribution of gravity data over the four main islands of Japan.

Digital elevation models (DEMs) covering the four main islands were prepared in 1999 by the GSI. The Hokkaido area is covered by a 250-m DEM, whereas the other three main islands (Honshu, Shikoku, and Kyushu) are covered by a 50-m DEM. Both the DEMs and gravity data coordinates are based on the Tokyo Datum, and so an initial datum transformation to the Japanese Geodetic Datum 2000 (JGD2000) was carried out for DEMs and the gravity data. Figure 3 shows the approximate topography of the study area.

Figure 3

Topography of the study area (units in m).

The method chosen for downward continuation of observed gravity on the topographical surface to the geoid plays a central role in precise geoid determination. Stokes's formula for the gravimetric geoid determination requires that no masses occur outside of the geoid and that the gravity anomaly should be referred to the geoid. This reduction was performed following the approach applied by Odera et al. (2012) to the previous geoid model for Japan. Although the classical Moritz formula (Moritz 1980) was used to estimate DTE in the previous geoid model for Japan, here, we evaluated DTE using the integral formula of Martinec and Vaníček (1994b). While both formulae use average topographic density, it should be noted that the classical Moritz formula is also based on the assumption that the distance between the computation and the running point is much greater than the height of the computation point (l _o  ≫ H _P ), which implies that it can only be used effectively for the far-zone integration area (Martinec et al. 1996). As a result, the effect of the near zone and the Bouguer shell, which cannot be derived in the planar model, is completely missing (e.g., Martinec and Vaníček 1994a; Nahavandchi 2000).

In an attempt to account for the contribution of the near zone, Martinec and Vaníček (1994b) divided the integration area (σ) into a near zone (σ₁) and a far zone (σ₂) based on a spherical approximation of the geoid. The resulting expression for DTE is given as (Martinec and Vaníček 1994b)

$$\begin{array}{l}\mathit{C}=-\frac{4\mathit{\pi }\mathit{G}{\mathit{\rho }}_0}{\mathit{R}}{\mathit{H}}_{\mathit{P}}^2+\frac{\mathit{G}{\mathit{\rho }}_0{\mathit{R}}^2}{2}{\displaystyle \underset{{\mathit{\sigma }}_1}{\iint }\frac{{\mathit{H}}^2-{\mathit{H}}_{\mathit{P}}^2}{{\mathit{l}}^3}\left(1-\frac{3{\mathit{H}}_{\mathit{P}}^2}{{\mathit{l}}^2}\right)}\mathit{d}\mathit{\sigma }\\ +\frac{\mathit{G}{\mathit{\rho }}_0{\mathit{R}}^2}{2}{\displaystyle \underset{{\mathit{\sigma }}_2}{\iint }\frac{{\mathit{H}}^2-{\mathit{H}}_{\mathit{P}}^2}{{\mathit{l}}_{\mathit{o}}^3}}\left(1-3{\text{sin}}^2\frac{\mathit{\psi }}{2}\right)\mathit{d}\mathit{\sigma },\end{array}$$

(3)

where R is the mean radius of the Earth, H _P is the height of the computation point, H is the height of the running point, σ is the surface integration element, ρ₀ is the average topographic density, and G is the Newtonian gravitational constant. The spatial distances l and l _o between the computation and the running points are given respectively as

$$\mathit{l}={\left({\mathit{r}}^2+{\mathit{R}}^2-2\mathit{rR}\text{cos}\mathit{\psi }\right)}^{\frac{1}{2}},$$

(4)

and

$${\mathit{l}}_{\mathit{o}}=2\mathit{R}\text{sin}\frac{\mathit{\psi }}{2},$$

(5)

where r is approximated as R + H _P .

The near zone is defined by the spherical distance (ψ₀) given by Martinec and Vaníček (1994b) as

$${\mathit{\psi }}_0=\sqrt{\frac{{\mathit{H}}_{\mathit{P}}}{\mathit{R}}}.$$

(6)

The corresponding PITE on the geoid (e.g., Martinec and Vaníček 1994a) is given as

$$\begin{array}{l}{\mathit{N}}_{\mathrm{\text{ind}}}=-\frac{2\mathit{\pi }\mathit{G}{\mathit{\rho }}_0{\mathit{H}}_{\mathit{P}}^2}{\mathit{\gamma }}+\\ +\frac{\mathit{G}{\mathit{\rho }}_0{\mathit{R}}^2}{\mathit{\gamma }}{\displaystyle \underset{\mathit{\sigma }}{\iint }\left[2\frac{\sqrt{{\mathit{l}}_{\mathit{o}}^2+{\mathit{H}}^2}-\sqrt{{\mathit{l}}_{\mathit{o}}^2+{\mathit{H}}_{\mathit{P}}^2}}{\mathit{R}}+\mathit{In}\frac{\frac{{\mathit{l}}_{\mathit{o}}}{2\mathit{R}}+\mathit{H}+\sqrt{{\mathit{l}}_{\mathit{o}}^2+{\mathit{H}}^2}}{\frac{{\mathit{l}}_{\mathit{o}}}{2\mathit{R}}+{\mathit{H}}_{\mathit{P}}+\sqrt{{\mathit{l}}_{\mathit{o}}^2+{\mathit{H}}_{\mathit{P}}^2}}-\frac{\mathit{H}-{\mathit{H}}_{\mathit{P}}}{{\mathit{l}}_{\mathit{o}}}\right]}\mathit{d}\mathit{\sigma }.\end{array}$$

(7)

The other applied reductions include a secondary indirect topographical effect (Vaníček et al. 1999), a second-order free-air reduction (Heiskanen and Moritz 1967), and an atmospheric correction (Wichiencharoen 1982a). The residual gravity anomalies are obtained from the reduced gravity anomalies on the geoid and gravity anomalies obtained from EGM2008 (degree 2,190 and order 2,159) in the three islands (Hokkaido, Honshu, and Kyushu), and GOCO02S/EGM2008 (degree 2,190 and order 2,159) in Shikoku. Kriging (Krige 1951) on a 1 × 1.5 arc-minute grid was performed for gridding of the residual gravity anomalies.

Geoid modeling over Japan

Although Stokes's integral formula for geoid determination (Stokes 1849; Heiskanen and Moritz 1967) was applied in this study, it requires that continuous gravity data be available and be used over the entire Earth. As this condition is not currently satisfied, a remove-compute-restore procedure was adopted here whereby a GGM provides the long-wavelength geoid undulations, a modified Stokes's formula provides the medium-wavelength components, and the indirect effect provides the short-wavelength components. In addition, the modified Stokes's kernel proposed by Meissl (1971) was used to minimize truncation errors. Hence, the expanded modified Stokes's formula, excluding small errors due to ellipsoidal effects, is given as

$$\mathit{N}={\mathit{N}}_{\mathrm{\text{GGM}}}+\frac{\mathit{R}}{4\mathit{\pi }\mathit{\gamma }}{\displaystyle \underset{{\mathit{\sigma }}_{\mathit{o}}}{\iint }\mathit{\Delta }{\mathit{g}}^{\mathit{r}}}{\mathit{S}}^{\mathrm{ME}}\left(\mathit{\psi }\right)\mathit{d}\mathit{\sigma }+{\mathit{N}}_{\mathrm{\text{ind}}},$$

(8)

where N is the gravimetric geoid undulation, N _GGM is the geoid undulation obtained from EGM2008 or GOCO02S/EGM2008 after applying the zero-degree term (with respect to GRS80), Δg^r is the residual gravity anomaly, N _ind is the indirect topographical effect, and S^ME(ψ) is Meissl's modified kernel.

The optimal spherical distance was evaluated empirically by comparing gravimetric and GPS/leveling geoid undulations at various spherical distances. These comparisons were performed at all 816 GPS/leveling points, and a standard deviation (S.D.) was computed for each spherical distance. As a result, a spherical cap size of 40 km was adopted for computations in this study because it gives the smallest S.D. when compared to the other spherical distances tested.

Results and discussion

Statistics relating to the topographical effects (computed by the integral formulae) and the corresponding values (computed by the classical formulae) are presented in Table 1. Computations were performed at the gravity data points.

Table 1 Statistics of the topographical effects over Japan

Figures 4 and 5 show the spatial distribution of DTE over Japan as computed by the classical Moritz formula and the integral formula of Martinec and Vaníček, respectively. Figures 6 and 7 show the spatial distribution of PITE over Japan as computed by a planar formula and the integral formula of Martinec and Vaníček, respectively. It should be noted that the spatial representation of the topographic effects over Japan was achieved by interpolation of the values from the gravity data points onto a 1 × 1.5 arc-minute grid for plotting, thus explaining why the minimum and maximum limits in Table 1 may vary slightly from those shown in Figures 4, 5, 6, 7.

Figure 4

DTE by classical Moritz formula over Japan (units in mGal).

Figure 5

DTE by integral formula of Martinec and Vaníček over Japan (units in mGal).

Figure 6

PITE by planar formula over Japan (units in cm).

Figure 7

PITE by integral formula of Martinec and Vaníček over Japan (units in cm).

The differences between gravimetric and GPS/leveling geoid undulations over Japan are shown in Figure 8, and the statistics of these comparisons are given in Table 2. The statistics of the differences between the previous gravimetric and GPS/leveling geoid undulations are also given in Table 2 in parentheses for comparison.

Figure 8

Differences between gravimetric and GPS/leveling geoid undulations over Japan (units in cm).

Table 2 Statistics of the differences between gravimetric and GPS/leveling geoid undulations over Japan and the six sub-regions

Table 2 shows that the West Honshu area has the smallest S.D. (±4.9 cm) between the gravimetric and GPS/leveling geoid undulations, while Central Honshu (a mountainous area) has the largest (±7.1 cm). When comparing these 816 GPS/leveling data points with the previous geoid model for Japan (Odera et al. 2012), there is a marked improvement in the overall S.D. from ±8.3 to ±7.5 cm.

It is evident that much of the improvement of the geoid is due to this study's treatment of the topography, given that the same data sets have been used for both models. For example, if EGM2008 is used over the four main islands (Hokkaido, Honshu, Shikoku, and Kyushu) and integral formulae are applied for the computation of topographic effects, the S.D. between gravimetric and GPS/leveling geoid undulations over the whole of Japan is ±7.6 cm, whereas it is ±8.7 cm over Shikoku only. Hence, the integral formulae of Martinec and Vaníček (1994a, b) seem to be more suitable for computing the DTE and PITE over Japan than the classical formulae. We also note the interesting effect of a deterioration of the S.D. over North Honshu as opportunity for further investigation. We suggest therefore that a significant improvement in the geoid model over Shikoku, from ±8.7 cm previously to ±6.6 cm currently, is mainly due to the contribution of GOCE data.

Conclusions

We described a procedure for the computation of an improved high-resolution gravimetric geoid model over Japan that exploits the currently available data sets. The integral formulae of Martinec and Vaníček (1994a, b) were used to compute the DTE and PITE over Japan for the first time and were found to perform better than the classical formulae over the studied area. A combination of GOCO02S and EGM2008 (GOCO02S/EGM2008) was used for geoid determination in the Shikoku area, while EGM2008 was used for geoid determination in the remaining islands (Hokkaido, Honshu, and Kyushu).

When comparing these 816 GPS/leveling data points with the previous gravimetric geoid model for Japan, a marked improvement in the overall S.D. from ±8.3 to ±7.5 cm was observed. It was noted that much of the improvement in the accuracy of the geoid over the study area was due to the use of the integral formulae for handling the topographical effects on gravity.

With the recent improvements in geoid determination techniques, the handling of the DTE and PITE needs further attention. However, the need for denser coverage and more accurate gravity data over Japan cannot be overemphasized. Despite this, the gravimetric geoid model obtained in this study performs significantly better than the previous versions and makes an important contribution towards the establishment of a geoid-consistent vertical datum over Japan.