SPHEROIDAL SPLINE INTERPOLATION AND ITS APPLICATION IN GEODESY

The aim of this paper is to study the spline interpolation problem in spheroidal geometry. We follow the minimization of the norm of the iterated Beltrami-Laplace and consecutive iterated Helmholtz operators for all functions belonging to an appropriate Hilbert space defined on the spheroid. By exploiting surface Green’s functions, reproducing kernels for discrete Dirichlet and Neumann conditions are constructed in the spheroidal geometry. According to a complete system of surface spheroidal harmonics, generalized Green’s functions are also defined. Based on the minimization problem and corresponding reproducing kernel, spline interpolant which minimizes the desired norm and satisfies the given discrete conditions is defined on the spheroidal surface. The application of the results in Geodesy is explained in the gravity data interpolation over the globe.


Introduction
This paper is concerned with finding a minimizer of the applied specific operator, which is either the iterated Beltrami-Laplace or consecutive iterated Helmholtz, to the particular Hilbert space in the sense of norm on a smooth manifold. The minimizer satisfies the Dirichlet or Neumann conditions at given points of the manifold.
The applications of this problem are in a number of fields in mathematical physics and engineering, e.g., Earth's gravity field and linear elasticity theory. The solution of the minimization problem has physical interpretation. For instance, in elasticity theory it means minimizing the bending energy of a thin layer.
In order to solve the minimization problem, spline interpolation approach is considered. Based on what the domain is and which boundary conditions are used, different problems of finding spline interpolant must be solved.
In the spherical domain, the spline interpolation has been investigated in a number of different studies including , Freeden (1981), Keller and Borkowski (2019), Sloan and Womersley (2002), Baramidze et al. (2006), and Wahba (1981, 1990. In , spherical splines and zonal kernels are defined in the context of Sobolev and reproducing kernel in interpolating the GRACE satellites' positions in their ground-track repeat orbit and the Total Electron Content (TEC) in the ionosphere are presented. In Klees et al. (2008), the concept of RBFs is employed to locally interpolate gravity data. Several types of spherical RBFs are introduced and implemented to interpolate residual gravity anomalies using GPS-Leveling data. In Sloan and Womersley (2002), the data interpolation in Geodesy is discussed by the concepts of minimum energy. Spherical harmonics are taken as the fundamental bases and the weight coefficients. Spherical interpolating and smoothing splines for a set of linearly independent evaluation functionals are discussed in Wahba (1981). Moreover, reproducing kernels are given in an analytic integral representation.
We are motivated by the minimization problem for data interpolation on the surface of a spheroid. Physical explanation for the importance of spheroidal surface can be observed in the interpolation of gravity data in Earth's gravity field, with its geometry being better defined with a spheroid. Earth is impacted by a number of different (inner and outer) forces, which make its shape irregular. To approximate this rough surface, we can use a sphere. However, observations have confirmed the spheroidal shape of the Earth. So, having a mathematical framework for interpolation on this special surface is highly needed. Although many works have been done on the spherical (spline) interpolation, there is little done on the spheroidal case. An important work for outer spheroidal spline, namely, Abel-Poisson kernel spline, has been done in Akhtar and Michel (2012); However, the case where data are on the surface of spheroid is fundamentally different and is investigated in the present paper.
This paper is organized as follows. Preliminaries and minimization problems are stated in Section 1. Section 2 is devoted to the surface Green's functions approach for Dirichlet and Neumann conditions at the given points of the spheroid. This section is ended with finding the iterated and consecutive iterated surface Green's functions. The definition of spline interpolant is developed in Section 3. Application of spheroidal spline interpolation for Gravity Data in Geodesy are presented in Section 4. At the end, conclusions are stated. In Appendix A, the eigenvalues and eigenfunctions of the Beltrami-Laplace equation are derived. Finally, the approximate formula for iterated Green's function is given in Appendix B.

Preliminaries and minimization problem
Similar to a sphere, an oblate spheroid is a surface of revolution with genus zero. Because of the inequality of its axes, it acquires a special geometry different from a sphere. Several important operators and definitions that are used frequently in the paper are given in the following.
Definition 1.1. Let a and b be the semi-major and semi-minor axes of the oblate spheroid, respectively. The eccentricity of the corresponding oblate spheroid is defined by (2) where, for , 1,2, , i j n =  , ij g and g are the elements and determinant of a metric tensor, respectively.
For the case 1 v = in (2) The i-th Helmholtz operator is defined as sum of the Beltrami-Laplace operator and the negative of its i-th eigenvalue i p , namely, The consecutive iterated Helmholtz operator of degree v , with its i-th element acting i q times, is defined as As an interpolation problem, the following set of J distinct points is given such that the given functionals applied to S at points 1 2 , , , J η η η  are known real values 1 2 , , , J U U U  , respectively. In this paper, the given functionals for the minimization problems (6) and (7) , are either of the following cases − discrete Dirichlet condition at , namely, where n is outward normal vector to  . Remark 1.2. The existence and uniqueness of the solutions of problems (6) and (7) are guaranteed (for more details, see Freeden, 1981 H  is of key importance to define the spline interpolant. This space can be written as the direct sum of the null space of the operator and the Hilbert space To find the spline interpolant in the Hilbert space ( ) = , has to be constructed.

Surface Green's functions
In order to define reproducing kernels, a common approach is to find surface Green's functions (see Freeden, 1981). According to the differential operators of the form (2)-(5), the surface Green's functions are derived for discrete Dirichlet and Neumann conditions. Following the method of Green's function, described in Greenberg (2015), it can be shown that

System of surface spheroidal harmonics for Green's functions
In order to derive Green's functions, one needs to have a fundamental system of surface spheroidal harmonics. As mentioned in Remark 1.3, we have the following decomposition If the solution has a zero degree eigenfunction, then the problem of finding the generalized Green's function has to be considered. The generalized Green's function is independent of the given data. Based on the definition of spline interpolant similar to Freeden (1984), Schreiner (2009), andWahba (1981), this type of Green's function can be exploited for both discrete Dirichlet and Neumann conditions. In the following subsections, the first Green's functions are introduced. These functions do not possess a zero order eigenfunction. With these functions, the fundamental system of surface harmonics are designed. Also the second Green's functions are created for discrete Neumann condition. It can be observed that the second Green's functions possess a zero order eigenfunction. In the case of generalized Green's functions, a set of complete orthonormal basis of surface spheroidal harmonics with a zero degree eigenfunction is formed.

First Green's function
For , ξ η∈ and  in the form of (2), the function ( ) , G ξ η is called the first Green's function or Green's function for Dirichlet condition, when it satisfies where δ denotes the Dirac delta function. Relation (10) holds for discrete Dirichlet condition (8). From (10) and (11), it is simply concluded that For finding the first Green's function, a general method, proposed in Felsen and Marcuritz (1994), is considered. According to this method, a singular point is placed on the spheroid's poles 0, θ = π . Thus the Green's function is singular at these points. Furthermore, the Green's function obeys the laws of wave propagation. For instance, the surface and its factor from which the wave crosses are indispensably important.
To solve (14), the method of separation of variables is considered as Based on the distribution theory (see details in Felsen and Marcuritz (1994)), the azimuth delta function can be written as By substituting (15) in (14) and (13), the first green's function is expressed by the following series From (16) and (13), it is deduced that Solving Equation (17) By integrating both sides of (17) in the interval and tending to zero, it is obtained that From Equation (19) and (18), The coefficient 0 a in (18) It can be shown that the Equation (20) leads to the eigenvalue problem The solutions of Equation (21) Therefore, m G can be written as Again by performing the same task as for (19), it can be concluded that  (22) and (23) make the series in (12) uniformly convergent.

Second Green's function
Again suppose that the operator  has the form (2) with 1 v = . The function ( ) , G ξ η is called second Green's function or Green's function for discrete Neumann condition, when it is the solution of the following problem is the surface area of the spheroid  . From (24) and (10), it is deduced that Remark 2.1. If the complete system of surface harmonics is chosen, then the generalized Green's function could be used for any set of linearly independent evaluation functionals (see Wahba, 1981). Therefore, the generalized Green's function under discrete Dirichlet conditions at  also satisfies (24).
In the following Proposition, the second and forth order approximate formulas are given for the second Green's function.
Proposition 2.2. The second and forth order expansions for the second Green's function of Beltrami-Laplace operator are, respectively, Proof. First of all, the following appropriate conformal mapping must be provided With this conformal map, the spheroid is mapped to the complex plane where two-dimensional variables can be represented in a complex variable. We set This means that the second Green's function is shown in the complex plane as the function ( ) Ø , ζ τ . By using (27) in (24) where ζ denotes the conjugate of ζ . By using (27) Choosing the positive root of Equation (32) Substituting (33) into (27) and (28)  The second green's function problem is closely associated with the Hamiltonian problem in spheroidal geometry. The Hamiltonian formula is a representation of the system's total dynamic energy for moving vortexes on the spheroid  . To achieve the Green's function, we change the Hamiltonian formula for dynamic energy (see more details in Castilho and Machado (2008)). For this reason, an appropriate conformal factor is proposed as Note that the conformal factor ( ) µ θ is not constant while the eccentricity satisfies 1 e < . Therefore, the function ( )  (2008), it is concluded that where j µ for 0,1, , j =  are coefficients of power series of function µ around e . In (35)  The surface area A  is smaller than the surface area of its corresponding sphere. Therefore, the spheroidal Green's function is smaller than the spherical Green's function.

Iterated Green's functions
We consider the iterated problem for the first and second Green's functions where the Beltrami-Laplace operator acts v times. From the general theory of Green's functions as presented in Freeden and Schreiner (2009), iterated Green's function v G satisfies the following convolution Similar to the proposed method in Greenberg (2015), this problem can be solved based on eigenvalue expansion. For solving (37) For the iterated problem with discrete Neumann condition at  , we get For The iterated problem with discrete Dirichlet conditions at  , we have

Green's function for the iterated Helmholtz operator
Suppose that  is of the form (4). The Green's function for the i-th Helmholtz operator is defined as By using (39) and (10), we have

Green's function for the consecutive iterated Helmholtz operator
We consider the consecutive Helmholtz operator by setting 0 1 i q q = = =  in (5). Based on the null space of the consecutive Helmholtz operator and following (39), the Green's function for the i-th degree satisfies Taking into account the following convolution and using the eigenvalue expansion procedure in Freeden and Schreiner (2009), it could be shown that (41) and (10), it is concluded that (38) and (40), the consecutive iterated Helmholtz Green's function reads as

Spline interpolant
So far, the Green's functions on the spheroid  have been derived. According to Freeden (1984) and (1981), we introduce the corresponding reproducing kernels and spline interpolants. For an admissible system of points { | 1, , } i i J η =  , we consider the Lagrange basis function where δ denotes the Kronecker symbol. Definition 3.1. For the unisolvent system of points and discrete Dirichlet and Neumann conditions -Reproducing kernel for the iterated Beltrami-Laplace operator is According to the considerations given in Freeden (1984), the definition of the spline function is presented in the following. , , where 1 J denotes the first 1 J elements of points which constitute an admissible system.
Note that the whole formulas for Green's functions and spline interpolants are valid and applicable in the prolate spheroidal coordinate system.

Application in gravity data interpolation in geodesy: interpolating global potential data
There are many applications to the spline functions, including gravity data interpolation, satellite position interpolation between observation points, and TEC interpolation between the discrete measurements in the Ionosphere (Keller & Borkowski, 2019). We show how the results obtained in the previous sections can be applied to globally interpolate gravity data. To make the 5° × 5° and 4° × 4° grid data, Potential values derived by EGM-2008 geopotential model (up to 2190 degree) are interpolated. Then these values are compared with the actual grid produced by the potential formula. The process of determining the unknown coefficients of the spline interpolation is exactly the same as that of the sphere (for more detail, see Freeden, 1981;Keller & Borkowski, 2019). We notice that the spherical coefficients in the potential formula have to convert to the ellipsoidal coefficients via the relations stated in Jekeli (1988). After the computations, we get the following Figure 1. Then the following steps are performed -According to (9), subtract the mean from the potential data on the surface of the ellipsoid. Consequently we get the following Figure 2. Figure 2. Potential over the globe, after subtracting the mean -By using the ellipsoidal spline interpolation, interpolate the data in the previous step for 4° × 4° and add the removed mean after interpolation. Then we get the following Figure 3. -To determine the difference, subtract the actual 4° × 4° grid derived from the potential formula from the 4° × 4° interpolated grid (Figure 4).

Conclusions
For solving the minimization problem on the spheroid, the surface Green's function approach has been chosen. Under the specific interpolation conditions, the Green's functions for the Beltrami-Laplace and consecutive Helmholtz operators and their iterations have been derived. By using the Green's functions, the corresponding reproducing kernels of the operators and the interpolation conditions on the given points have been introduced. Based on the reproducing kernels, the spline interpolants, as the minimizer of the problem, are presented. It is verified that this work is the generalization of the minimization problem on a sphere. This means that the spline interpolant formulas on a sphere could be obtained by the corresponding formulas on a spheroid when the eccentricity tends to zero. This work can be used in different areas of study, including Earth's gravity field where the geometrical structure of the Earth is better modeled with a spheroid. In future work, we intend to extend the minimization problem of spheroidal smoothing spline to ellipsoidal geometry for a set of linearly independent evaluation functionals.

Author contributions
All four authors contributed to the derivation of the formulae and the writing of the manuscript. A. Safari and M. Kiani did the data analysis for the section in numerical computations.
Based on Lagrange solution of an inhomogeneous second order differential equation, the solution of (47), expressed in an implicit form, reads as Equation (50) with 0 e = is the Legendre equation. Based on the series expansion method, a Taylor series representation around zero exists for Equation (50). Therefore, the solution of Equation (50)