An Online Gravity Modeling Method Applied for High Precision Free-INS

For real-time solution of inertial navigation system (INS), the high-degree spherical harmonic gravity model (SHM) is not applicable because of its time and space complexity, in which traditional normal gravity model (NGM) has been the dominant technique for gravity compensation. In this paper, a two-dimensional second-order polynomial model is derived from SHM according to the approximate linear characteristic of regional disturbing potential. Firstly, deflections of vertical (DOVs) on dense grids are calculated with SHM in an external computer. And then, the polynomial coefficients are obtained using these DOVs. To achieve global navigation, the coefficients and applicable region of polynomial model are both updated synchronously in above computer. Compared with high-degree SHM, the polynomial model takes less storage and computational time at the expense of minor precision. Meanwhile, the model is more accurate than NGM. Finally, numerical test and INS experiment show that the proposed method outperforms traditional gravity models applied for high precision free-INS.


Introduction
Inertial navigation systems (INSs) are employed throughout all branches of the military and in many civil platforms, due to their significant advantages of providing continuous position, velocity and attitude information, and being invulnerable to external interference [1]. In the procedure of real-time navigation, the gravity acceleration is one of the significant pieces of external reference information needed for INSs [2]. This continuous gravity information is routinely obtained from the normal gravity model (NGM). With the development of inertial sensors, the availability of high-precision gravity information has become a crucial issue for preparing accurate INSs [3]. Therefore, a more accurate gravity description in high-precision INSs should be promoted. Besides online gravity models, en route gravity compensation methods using gradiometers are reported [4,5]. However, due to its independence from the accuracy of gravity sensitive instrumentation, gravity modeling has always been introduced into INSs [5,6].
As mentioned above, the normal gravity models (NGMs) described on the surface of reference ellipsoids are conventionally used to compensate for gravitational effects on accelerometers [7]. These models are usually expressed as simple combinations of first-and second-order sine functions. They are successfully applied for current low or medium accuracy INSs. However, for high-precision systems, the NGMs are inefficient. To increase NGMs' accuracy, researchers have developed several improved NGMs [8][9][10], while the residual error of the normal gravity vector to the actual one (i.e., gravity disturbance vector, GDV) is still ignored. In the field of geodesy, spherical harmonic series, commonly used to represent the Earth's gravitational field, are routinely expanded to ultra-high degree (>2160). In addition, the spherical harmonic gravity model (SHM) could describe the GDV accurately. That promotes the possibility of the application of spherical harmonic gravity model (SHM) for INSs to compensate for GDV-induced errors. However, we have demonstrated that only SHM of degree 12 or less is suitable for high frequency real-time solution of INSs, which still retains on average a 27.84 mGal horizontal gravity error globally [11]. Since SHM's computational time increases quadratically with its degree, it greatly exceeds a single navigation period. According to our surveys, no studies have mentioned the application of high-degree SHM for INSs. Even so, the high-fidelity SHM has always attracted the attentions of scientists in astrodynamics to develop efficient models for rapid and precise orbit propagation.
To improve the efficiency of SHMs' execution, 3D interpolation models perform more acceptably for the Earth's application [12]. The modeling method, effectively a trade-off between speed and memory, stores nodal gravity field information on an equivalent spherical grid. This information is previously calculated with SHM on a computer. Then, a variety of interpolation methods are employed to solve the gravity for real-time orbit determination. Junkins first proposed this modeling approach in 1976. However, it was limited by the low quality of space-borne computers at the time, because it needs to store almost 30,000 coefficients [13]. In 1996, Hujsak introduced the concept of pseudocenter and achieved a 100-fold speedup but with consumption of precision compared to SHM of degree 70. The coordinates of the pseudocenter can be used to derive gravity information, which reduces the memory requirement to the order of several megabytes [14]. Recently, the developments of memory and processor technology applied for satellites have motivated more works on the 3D interpolation models. Interpolation methods, such as weighting functions, wavelets, B-splines and octrees, were employed in the models to reduce the computational complexity [12,15]. Efficient models, such as the cubed-sphere model and the fetch model, were established to balance accuracy with complexity. However, in principle these 3D models improve the calculation speed at the expense of storage, which is the most limited resource for navigation computers, contrary to space satellites. Therefore, for high-precision INS, the direct application of these efficient models will be limited.
With the motivation of compensating the GDV-induced error for accurate INSs, the high-fidelity SHM is an obvious choice. To apply the model for real-time navigation, especially for navigation in rugged regions, we will modify it as an online gravity model in this paper. Therefore, considering the advantages of 3D interpolation models, we propose a more time-and space-efficient online modeling method. The method combines a regional two-dimensional second-order polynomial model (PM) with an online strategy for the coefficients. The PM is selected to establish gravity on a 6 × 6 navigation area, where the components of GDV vary almost linearly with latitude and longitude. To achieve global navigation, the coefficients of the PM update as the vehicle gets close to the boundary of the area. The updating of coefficients and the switch of navigation area are executed in an external gravity computer, where the PM is fitted to nodal gravity information calculated by SHM. This design is superior to data-based compensation methods that occupy several times the memory resources (see [8]) of harmonic coefficients. Compared with 3D interpolation models, the proposed model possesses only twelve coefficients that require a few bytes of memory. Moreover, it performs more efficiently relative to the interpolation methods.
The outline of this paper is as follows: Section 1 is the Introduction; the definition of gravity disturbance vector and its expression as spherical harmonic series are reviewed briefly in Section 2; Section 3 gives a demonstration of variations of gravity components on a small-scale range, and a numerical test is carried out to verify the theoretical derivation. In Section 4, the strategy of the novel online modeling method and gravity compensation is designed. Experimental results and discussions of the proposed method applied for a real free-INS are shown in Section 5, and conclusions are presented in Section 6.

Definition of Gravity Disturbance Vector
As Figure 1 shows, a point P on the geoid is projected onto the point Q on the surface of the ellipsoid by means of the ellipsoidal normal. The gravity disturbance vector (GDV) is defined as the difference between the gravity vector g at P and the normal gravity vector γ at the same point P [7]: The difference in magnitude is the gravity disturbance, which historically presents as gravity anomaly ∆g, for it is more available and being processed than gravity disturbance:

Definition of Gravity Disturbance Vector
As Figure 1 shows, a point P on the geoid is projected onto the point Q on the surface of the ellipsoid by means of the ellipsoidal normal. The gravity disturbance vector (GDV) is defined as the difference between the gravity vector g at P and the normal gravity vector γ at the same point P [7]: The difference in magnitude is the gravity disturbance, which historically presents as gravity anomaly Δg, for it is more available and being processed than gravity disturbance: The difference in direction is the deflection of vertical (DOV), which has two components, a north-south component ξ and an east-west component η. Therefore, the GDV on the Earth's surface for a high precision INS can be expressed as: where, the superscript n denotes gravity defined in local geodetic frame. γ is the magnitude of the normal gravity vector γ, whose horizontal components on the ellipsoid's surface are zero, i.e., γ = [0 0 -γ]. The GDV and normal gravity vector constitute the gravity vector compensated in high precision INS. The expressions of GDV's three components obviously determine the accuracy and efficiency of the gravity model.

Spherical Harmonic Model
Commonly, the spherical harmonic expansion is used to represent the Earth's gravitational field with high fidelity. The Earth Gravitational Model 2008 (EGM2008) computed by the U.S. National Geospatial-Intelligence Agency (NGA) is complete to degree 2160 [16]. It allows computing the disturbing potential and a number of gravity field related quantities, such as height anomalies, gravity anomalies and vertical deflections to a spatial resolution of about five arc minutes. Therefore, it is an optimal choice to invite the SHM for describing the GDV. Firstly, we write the expression of disturbing potential, of which the GDV is the gradient vector: where, f is Newton's gravitational constant; M is mass of the earth; ρ is the geocentric radius distance; θ is spherical polar angle; λ is the geocentric longitude; C′nm and Snm are constant coefficients of the spherical harmonic of degree n and order m, among them C′nm is modified by reference ellipsoid parameters; Pnm is the normalized associated Legendre function of the first kind. Then, we have the components of GDV on a computation point (situated at the Earth's surface) as follows: The difference in direction is the deflection of vertical (DOV), which has two components, a north-south component ξ and an east-west component η. Therefore, the GDV on the Earth's surface for a high precision INS can be expressed as: where, the superscript n denotes gravity defined in local geodetic frame. γ is the magnitude of the normal gravity vector γ, whose horizontal components on the ellipsoid's surface are zero, i.e., γ = [0 0 -γ]. The GDV and normal gravity vector constitute the gravity vector compensated in high precision INS. The expressions of GDV's three components obviously determine the accuracy and efficiency of the gravity model.

Spherical Harmonic Model
Commonly, the spherical harmonic expansion is used to represent the Earth's gravitational field with high fidelity. The Earth Gravitational Model 2008 (EGM2008) computed by the U.S. National Geospatial-Intelligence Agency (NGA) is complete to degree 2160 [16]. It allows computing the disturbing potential and a number of gravity field related quantities, such as height anomalies, gravity anomalies and vertical deflections to a spatial resolution of about five arc minutes. Therefore, it is an optimal choice to invite the SHM for describing the GDV. Firstly, we write the expression of disturbing potential, of which the GDV is the gradient vector: where, f is Newton's gravitational constant; M is mass of the earth; ρ is the geocentric radius distance; θ is spherical polar angle; λ is the geocentric longitude; C nm and S nm are constant coefficients of the spherical harmonic of degree n and order m, among them C nm is modified by reference ellipsoid parameters; P nm is the normalized associated Legendre function of the first kind. Then, we have the components of GDV on a computation point (situated at the Earth's surface) as follows: Accordingly, all components of the gravity model are expressed in the form of spherical harmonic series. This will impose a large computational complexity on the navigation solution. To retain the advantage of high-degree SHM and use it for a high precision INS, a simplified expression with less demands for accuracy should be investigated. Our following work will concentrate on these issues.

Gravity Disturbing Potential in Small Area
Since gravity disturbing potential formulation can be used to derive variations in gravity anomaly and vertical deflections, our discussions focus on Equation (4). When truncated at degree n max , the spatial resolution of the spherical harmonic expansion can be calculated as follows: In the case of EGM2008, the maximum degree n max of 2160 implies a resolution of 5 arc minutes. Gravity disturbing potential of degree n is expressed as: The quantity of value of zero is at most n and 2n in the directions of θ and λ, respectively. The interval of these zeros in both directions is equal to 180 • /n, i.e., the spatial resolution. That means the highest frequency of T n is n/2π. Figure 2 shows one period of function P n (cosθ) in the direction of θ.
Accordingly, all components of the gravity model are expressed in the form of spherical harmonic series. This will impose a large computational complexity on the navigation solution. To retain the advantage of high-degree SHM and use it for a high precision INS, a simplified expression with less demands for accuracy should be investigated. Our following work will concentrate on these issues.

Gravity Disturbing Potential in Small Area
Since gravity disturbing potential formulation can be used to derive variations in gravity anomaly and vertical deflections, our discussions focus on Equation (4). When truncated at degree nmax, the spatial resolution of the spherical harmonic expansion can be calculated as follows: In the case of EGM2008, the maximum degree nmax of 2160 implies a resolution of 5 arc minutes. Gravity disturbing potential of degree n is expressed as: The quantity of value of zero is at most n and 2n in the directions of θ and λ, respectively. The interval of these zeros in both directions is equal to 180°/n, i.e., the spatial resolution. That means the highest frequency of Tn is n/2π. Figure 2 shows one period of function Pn(cosθ) in the direction of θ. In any 180°/n × 180°/n region of the sphere, the average gravity field is usually used for spherical harmonic synthesis instead of discrete values [17][18][19]. Therefore, we take the sense that some regular features would be explored from the gravity calculated with SHM in small size area. Here, we re-write the gravity disturbing potential as [20]: with the lumped coefficients: In any 180 • /n × 180 • /n region of the sphere, the average gravity field is usually used for spherical harmonic synthesis instead of discrete values [17][18][19]. Therefore, we take the sense that some regular features would be explored from the gravity calculated with SHM in small size area. Here, we re-write the gravity disturbing potential as [20]: with the lumped coefficients: a ρ n C nm P nm (cosθ) a ρ n S nm P nm (cosθ) (11) In spherical harmonic analysis and synthesis, the surface of the Earth is gridded as numbers of subblock ∆θ × ∆λ [19]. We describe the disturbing potential values on a global grid as follows: Then, we have its neighbor disturbing potential on the same latitude: Considering the variation between the two points, we make a subtraction: 2 c m cos m l + 1 2 ∆λ + s m sin m l + 1 2 ∆λ sin m∆λ 2 (14) This difference oscillates at the frequency m/4π, which will vary with the quantity of m. The highest frequency is at the point m = 2160, in the case of EGM2008. Thus, the level of nonlinearity for this difference is determined by the sum of sin(m∆λ/2), where m ∈ [0, 2160]. Then, we turn our focus to the grid interval ∆λ. If it is less than the model's resolution, i.e., ∆λ < ∆Ψ, the value of m∆λ/2 will be less than π/2 at m = 2160. That means the sine function will be a small value that could be simply treated as a linear or low order polynomial, when ∆λ = ∆Ψ/2. For the sum of sin(m∆λ/2), all of the calculated angles are less than π/4. Assuming the coefficients of sine functions to be a series of normalized random data in the range of [−1, 1], the sum of sin(m∆λ/2) is shown in Figure 3 coupled with its fitted polynomial functions of different order. The inset of Figure 3 shows the differences between the regressed functions and the basic function.
In spherical harmonic analysis and synthesis, the surface of the Earth is gridded as numbers of subblock Δθ × Δλ [19]. We describe the disturbing potential values on a global grid as follows: Then, we have its neighbor disturbing potential on the same latitude: Considering the variation between the two points, we make a subtraction: This difference oscillates at the frequency m/4π, which will vary with the quantity of m. The highest frequency is at the point m = 2160, in the case of EGM2008. Thus, the level of nonlinearity for this difference is determined by the sum of sin(mΔλ/2), where m  [0, 2160]. Then, we turn our focus to the grid interval Δλ. If it is less than the model's resolution, i.e., Δλ < ΔΨ, the value of mΔλ/2 will be less than π/2 at m = 2160. That means the sine function will be a small value that could be simply treated as a linear or low order polynomial, when Δλ = ΔΨ/2. For the sum of sin(mΔλ/2), all of the calculated angles are less than π/4. Assuming the coefficients of sine functions to be a series of normalized random data in the range of [−1, 1], the sum of sin(mΔλ/2) is shown in Figure 3 coupled with its fitted polynomial functions of different order. The inset of Figure 3 shows the differences between the regressed functions and the basic function. It is observed that the residual error of a third-order polynomial function fitted to the sum of sine functions is nearly 0.02. This is less than 0.2 percent of the peak value. The polynomial functions of order greater than three performed better when the interval Δλ is in the range of [−2.5, 2.5] arc minutes. Similarly, we could obtain the same results on the direction of Δθ according to the property of normalized associated Legendre function Pnm(cosθ). Therefore, we propose that two-dimensional polynomials (for the model is established on the Earth's surface, we consider the coordinate in the vertical direction to be constant here) of low order could be used to fit to any components of GDV in the area with size of half of the SHM's resolution. Numerical tests address this in the following section. It is observed that the residual error of a third-order polynomial function fitted to the sum of sine functions is nearly 0.02. This is less than 0.2 percent of the peak value. The polynomial functions of order greater than three performed better when the interval ∆λ is in the range of [−2.5, 2.5] arc minutes. Similarly, we could obtain the same results on the direction of ∆θ according to the property of normalized associated Legendre function P nm (cosθ). Therefore, we propose that two-dimensional polynomials (for the model is established on the Earth's surface, we consider the coordinate in the vertical direction to be constant here) of low order could be used to fit to any components of GDV in the area with size of half of the SHM's resolution. Numerical tests address this in the following section.

Test Design
To illustrate the theory above, two test areas are chosen, one in the North China Plain and another in the Himalayas. Based on EGM2008 and maximum degree 2160, 900 and 360, basic values of DOVs are generated on dense grids with 50 × 50 points. In addition, the database is obtained on three subregions of each test area, sizes of which are corresponding to 15 × 15 , 6 × 6 and 2.5 × 2.5 . A remark made here that we discuss the horizontal components of GDV only, because the vertical channel of INS normally requires some external aiding, such as barometric altitude [21].
The 50 × 50 DOVs of each area are shown in Tables 1 and 2. Recollecting DOVs on 15 × 15 knots of above areas, we fit two-dimensional second-and third-order polynomials to them. And then, DOVs on each 50 × 50 grids generated from PMs could be used to compare with the above database.

Test Results
Subtracting DOVs calculated by the PMs from database, we get the root mean square (RMS), the maximum and minimum values of the approximation errors. Table 3 reports the descriptive  statistics for the North China Plain and Table 4 for the Himalaya region.
In both test areas, the approximation errors show the same level of accuracy for subregions whose sizes are equal to half of the SHM's resolution. Specifically, using a second-order polynomial produces the RMS-errors of 10 −2 arc·s in smooth area and 10 −1 arc·s in rugged area. A third-order polynomial reduces the approximation errors to 10 −3 and 10 −2 arc·s, respectively. As the size is reduced, the approximation accuracy improves.
Actually, a 0.1 arc·s DOV induces a position error on the order of several meters per hour [22]. Therefore, to involve a larger applicative area of the fitted polynomial, we set this 0.1 arc·s DOV as a permission index of model approximation. Further, this index could theoretically satisfy the accuracy demand of a modern high-precision navigation system.

Online Gravity Modeling in INS
Since the two-dimensional polynomial of low order represents the regional gravity precisely, we propose the online gravity modeling scheme shown in Figure 4.

Test Results
Subtracting DOVs calculated by the PMs from database, we get the root mean square (RMS), the maximum and minimum values of the approximation errors. Table 3 reports the descriptive  statistics for the North China Plain and Table 4 for the Himalaya region.
In both test areas, the approximation errors show the same level of accuracy for subregions whose sizes are equal to half of the SHM's resolution. Specifically, using a second-order polynomial produces the RMS-errors of 10 −2 arc·s in smooth area and 10 −1 arc·s in rugged area. A third-order polynomial reduces the approximation errors to 10 −3 and 10 −2 arc·s, respectively. As the size is reduced, the approximation accuracy improves.
Actually, a 0.1 arc·s DOV induces a position error on the order of several meters per hour [22]. Therefore, to involve a larger applicative area of the fitted polynomial, we set this 0.1 arc·s DOV as a permission index of model approximation. Further, this index could theoretically satisfy the accuracy demand of a modern high-precision navigation system.

Online Gravity Modeling in INS
Since the two-dimensional polynomial of low order represents the regional gravity precisely, we propose the online gravity modeling scheme shown in Figure 4. As it shown, our navigation system comprises the navigation computer and the gravity computer [2,23]. The online modeling method provides PM for dynamic INS in real-time and global navigation. To update PM for varying navigation areas, the method will be aided by the external gravity computer. Once the vehicle is close to the boundary of the current area, the coefficients of PM and next applicable area will be updated in the external computer. The modeling scheme is implemented as follows: Step 1: INS gives the motion information of vehicle to the external computer, consisting of position, velocity and acceleration.
Step 2: The prediction module calculates the applicative area of the updated model, concluding maximum latitude, longitude and minimum latitude, longitude. The prime parameters for model configuration in the external computer are the degree of SHM, size of area and grid density, which we will discuss in the following section. As it shown, our navigation system comprises the navigation computer and the gravity computer [2,23]. The online modeling method provides PM for dynamic INS in real-time and global navigation. To update PM for varying navigation areas, the method will be aided by the external gravity computer. Once the vehicle is close to the boundary of the current area, the coefficients of PM and next applicable area will be updated in the external computer. The modeling scheme is implemented as follows: Step 1: INS gives the motion information of vehicle to the external computer, consisting of position, velocity and acceleration.
Step 2: The prediction module calculates the applicative area of the updated model, concluding maximum latitude, longitude and minimum latitude, longitude. The prime parameters for model configuration in the external computer are the degree of SHM, size of area and grid density, which we will discuss in the following section.

Degree Selection of Base Model
The SHM as the base model has two types of errors: commission errors and omission errors. The commission errors are due to incorrectly estimating coefficients of the model, and the omission errors are due to truncation errors of the short wavelength variations of the gravity potential which are not modeled [15,24]. The degree variances of DOVs and error degree variances are used to describe these errors, whose expressions are as follows: In the case of EGM2008, the DOVs degree variances, error degree variances and their ratio are plotted in Figure 5.

Degree Selection of Base Model
The SHM as the base model has two types of errors: commission errors and omission errors. The commission errors are due to incorrectly estimating coefficients of the model, and the omission errors are due to truncation errors of the short wavelength variations of the gravity potential which are not modeled [15,24]. The degree variances of DOVs and error degree variances are used to describe these errors, whose expressions are as follows:     In the case of EGM2008, the DOVs degree variances, error degree variances and their ratio are plotted in Figure 5. According to Figure 5, the signal to noise ratio (SNR) becomes 1 near degree 1840 and 0.1 near degree 900. Compared with SHM of degree 2160, the commission errors and omission errors in Figure 6 illustrate the smallest value in degree 820. Combining with our previous research that only one-third omission errors occupied by degree from 141 to 2160 [10], we choose SHM of degree 820 as our base model here. This model has a SNR of 0.08, and will guarantee the computational efficiency at the expense of a little accuracy. In this case, the commission error is 0.53 arc·s. In addition with the omission error of 0.08 arc·s, it yields a total error of 0.61 arc·s. Therefore, added to the polynomial approximation error of 0.1 arc·s, the 0.71 arc·s error will induce position error ranging from ten meters to thirty meters [22]. For an INS with precision demand 0.2 nm/h, the PM-induced residual error accounts for less than ten percent of the total position error.

Configuration of Polynomial Model
As mentioned in the preceding section, the main focus of our modeling method is to reconfigure SHM of degree 820 with a low order polynomial for local areas. The resolution of the model is about 13.2 arc·min. According to the analysis of Section 3, the size of the local area should not more than 6.6 arc min. To facilitate the configuration, we round this value down to 6 arc·min.
Subdividing the 6′ × 6′ area into several subintervals with equal angle along latitude and longitude, we obtain DOVs on the grid calculated by SHM. These basic values are utilized to configure our polynomial model.  According to Figure 5, the signal to noise ratio (SNR) becomes 1 near degree 1840 and 0.1 near degree 900. Compared with SHM of degree 2160, the commission errors and omission errors in Figure 6 illustrate the smallest value in degree 820. Combining with our previous research that only one-third omission errors occupied by degree from 141 to 2160 [10], we choose SHM of degree 820 as our base model here. This model has a SNR of 0.08, and will guarantee the computational efficiency at the expense of a little accuracy. In this case, the commission error is 0.53 arc·s. In addition with the omission error of 0.08 arc·s, it yields a total error of 0.61 arc·s. Therefore, added to the polynomial approximation error of 0.1 arc·s, the 0.71 arc·s error will induce position error ranging from ten meters to thirty meters [22]. For an INS with precision demand 0.2 nm/h, the PM-induced residual error accounts for less than ten percent of the total position error.

Configuration of Polynomial Model
As mentioned in the preceding section, the main focus of our modeling method is to reconfigure SHM of degree 820 with a low order polynomial for local areas. The resolution of the model is about 13.2 arc·min. According to the analysis of Section 3, the size of the local area should not more than 6.6 arc·min. To facilitate the configuration, we round this value down to 6 arc·min.
Subdividing the 6 × 6 area into several subintervals with equal angle along latitude and longitude, we obtain DOVs on the grid calculated by SHM. These basic values are utilized to configure our polynomial model. On the other hand, except for accuracy, the execution time of model configuration dominates our designation. A proper order of polynomial and knot density of area should be analyzed in modeling procedure. They are the primary limitations of the execution time for model preparation.
To investigate these values, we establish PM on the subdivision as Figure 7 shows: The basic data chosen here is DOVs from two 6′ × 6′ areas in Section 3.2, which are calculated by the base model. We fit first-to fifth-order PMs to the grid DOVs, whose number is k × k. In the test, the value of k is set as 3, 6, 9 and 12 respectively. The performances of PMs are tested through results of ten million points in each test area. The descriptive statistics (lg(RMS)) of the polynomial approximate results minus basic data are reported in Figure 8. On the other hand, except for accuracy, the execution time of model configuration dominates our designation. A proper order of polynomial and knot density of area should be analyzed in modeling procedure. They are the primary limitations of the execution time for model preparation.
To investigate these values, we establish PM on the subdivision as Figure 7 shows: The basic data chosen here is DOVs from two 6′ × 6′ areas in Section 3.2, which are calculated by the base model. We fit first-to fifth-order PMs to the grid DOVs, whose number is k × k. In the test, the value of k is set as 3, 6, 9 and 12 respectively. The performances of PMs are tested through results of ten million points in each test area. The descriptive statistics (lg(RMS)) of the polynomial approximate results minus basic data are reported in Figure 8. The basic data chosen here is DOVs from two 6 × 6 areas in Section 3.2, which are calculated by the base model. We fit first-to fifth-order PMs to the grid DOVs, whose number is k × k. In the test, the value of k is set as 3, 6, 9 and 12 respectively. The performances of PMs are tested through results of ten million points in each test area. The descriptive statistics (lg(RMS)) of the polynomial approximate results minus basic data are reported in Figure 8. On the other hand, except for accuracy, the execution time of model configuration dominates our designation. A proper order of polynomial and knot density of area should be analyzed in modeling procedure. They are the primary limitations of the execution time for model preparation.
To investigate these values, we establish PM on the subdivision as Figure 7 shows: The basic data chosen here is DOVs from two 6′ × 6′ areas in Section 3.2, which are calculated by the base model. We fit first-to fifth-order PMs to the grid DOVs, whose number is k × k. In the test, the value of k is set as 3, 6, 9 and 12 respectively. The performances of PMs are tested through results of ten million points in each test area. The descriptive statistics (lg(RMS)) of the polynomial approximate results minus basic data are reported in Figure 8.  In Figure 8, the contour label-2 denotes the RMS of errors approximated with PM is 0.01 arc·s. According to the result of fitting a second-order polynomial to the 3 × 3 grid DOVs, the RMS error could reduce to less than 0.1 arc·s in both the smooth area and rugged area. Therefore, when we choose the order of polynomial as 2 and the grid density as 3, the accuracy requirement of INS could be fulfilled greatly.
Finally, our approximation formulae of DOVs expressed as second-order polynomial of two dimensions are as follows: ξ assumes the form: η assumes the form: The coefficients aij, and bij (i, j = 0,1,2) will be determined by disturbance components at 3 × 3 points of predicted 6′ × 6′ navigation area. An additional topic that has to be presented here is the balance between size of area and grid density. During the update of model coefficients in the external computer, we hope it takes as little computational time as possible. Accordingly, we should choose a small value of k, based on which the size of area should be as large as possible. Taking the two test areas for example, we calculated the RMS results of errors for different size of area, when k = 3 and 6. As Figure 9 shows, a 6′ × 6′ update area of model is the best choice for the approximate accuracy of 0.1 arc·s.  In Figure 8, the contour label-2 denotes the RMS of errors approximated with PM is 0.01 arc·s. According to the result of fitting a second-order polynomial to the 3 × 3 grid DOVs, the RMS error could reduce to less than 0.1 arc·s in both the smooth area and rugged area. Therefore, when we choose the order of polynomial as 2 and the grid density as 3, the accuracy requirement of INS could be fulfilled greatly.
Finally, our approximation formulae of DOVs expressed as second-order polynomial of two dimensions are as follows: ξ assumes the form: η assumes the form: The coefficients a ij , and b ij (i, j = 0,1,2) will be determined by disturbance components at 3 × 3 points of predicted 6 × 6 navigation area. An additional topic that has to be presented here is the balance between size of area and grid density. During the update of model coefficients in the external computer, we hope it takes as little computational time as possible. Accordingly, we should choose a small value of k, based on which the size of area should be as large as possible. Taking the two test areas for example, we calculated the RMS results of errors for different size of area, when k = 3 and 6. As Figure 9 shows, a 6 × 6 update area of model is the best choice for the approximate accuracy of 0.1 arc·s.  In Figure 8, the contour label-2 denotes the RMS of errors approximated with PM is 0.01 arc·s. According to the result of fitting a second-order polynomial to the 3 × 3 grid DOVs, the RMS error could reduce to less than 0.1 arc·s in both the smooth area and rugged area. Therefore, when we choose the order of polynomial as 2 and the grid density as 3, the accuracy requirement of INS could be fulfilled greatly.
Finally, our approximation formulae of DOVs expressed as second-order polynomial of two dimensions are as follows: ξ assumes the form: η assumes the form: The coefficients aij, and bij (i, j = 0,1,2) will be determined by disturbance components at 3 × 3 points of predicted 6′ × 6′ navigation area. An additional topic that has to be presented here is the balance between size of area and grid density. During the update of model coefficients in the external computer, we hope it takes as little computational time as possible. Accordingly, we should choose a small value of k, based on which the size of area should be as large as possible. Taking the two test areas for example, we calculated the RMS results of errors for different size of area, when k = 3 and 6. As Figure 9 shows, a 6′ × 6′ update area of model is the best choice for the approximate accuracy of 0.1 arc·s.

Computational Complexity
The computational complexity of a model itself should be of concern except for its precision. In this section, we analyze the time and space complexity of the proposed modeling method. The arithmetic operations involved include addition (ADD), multiplication (MUL), division (DIV) and comparison (COMP). Considering the procedure of model configuration mentioned above, our online modeling method concludes two parts: one is the real-time computation of gravity vector for free-INS, and the other is the synchronous update of PM's coefficients in external computer.
It is obvious that there is two second-order polynomials and twelve coefficients in the real-time model computation. Therefore, the total operations of real-time model computation are 10ADD+18MUL, and the space complexity is 12 float. Storage taken by the model in navigation computer is approximate 0.09 KB. Tests for timing and accuracy reported in Table 5 were done by comparing performance of PM to routine NGM, modified NGM proposed in reference [8], the Cubed-Sphere model and SHM of degree 12 and 820. The test consists of computing gravity vectors for 10,000 randomly generated points, then comparing average execution times and accuracy. The speed-up factor in the fifth column is obtained by dividing the execution time for the SHM of degree 820 by the execution times for other models. The l2-norm of the difference between the gravity vectors produced by the first five models and those of SHM of degree 820 are reported, whose largest errors are in the sixth column. Observe that, the proposed PM performs outstandingly in both accuracy and execution time.
For the operation of updating the model's coefficients, computations of nine points with SHM of degree 820 are needed. Because the cm, sm are independent of λ, Equation (11) needs to be evaluated only once for computation points densely spaced along a parallel (i.e., φ constant), and the efficiency can be further increased for points equally spaced in longitude [17]. According to the time complexity of a single point analyzed in our previous work, the operations of the 3 × 3 grids are approximate to (13ADD+18MUL+3DIV+1COMP)(nmax + 2)(nmax − 1), where nmax is equal to 820. Compared with this process of gravity computation, the execution time of the fit algorithm to discrete data could be ignored. Without additional stored parameters, the space complexity of execution of coefficients update is equal to 2(nmax + 4)(nmax − 1) float. When nmax = 820, the storage occupied is 11.4 MB. Toward a processor with 2.6 GHz master frequency, the computation time is on average approximately 2.54 s. It means that the prediction of navigation area should be finished to retain the update time for model coefficients before the vehicle reaches the boundary of the present navigation area.

Computational Complexity
The computational complexity of a model itself should be of concern except for its precision. In this section, we analyze the time and space complexity of the proposed modeling method. The arithmetic operations involved include addition (ADD), multiplication (MUL), division (DIV) and comparison (COMP). Considering the procedure of model configuration mentioned above, our online modeling method concludes two parts: one is the real-time computation of gravity vector for free-INS, and the other is the synchronous update of PM's coefficients in external computer.
It is obvious that there is two second-order polynomials and twelve coefficients in the real-time model computation. Therefore, the total operations of real-time model computation are 10ADD+18MUL, and the space complexity is 12 float. Storage taken by the model in navigation computer is approximate 0.09 KB. Tests for timing and accuracy reported in Table 5 were done by comparing performance of PM to routine NGM, modified NGM proposed in reference [8], the Cubed-Sphere model and SHM of degree 12 and 820. The test consists of computing gravity vectors for 10,000 randomly generated points, then comparing average execution times and accuracy. The speed-up factor in the fifth column is obtained by dividing the execution time for the SHM of degree 820 by the execution times for other models. The l2-norm of the difference between the gravity vectors produced by the first five models and those of SHM of degree 820 are reported, whose largest errors are in the sixth column. Observe that, the proposed PM performs outstandingly in both accuracy and execution time.
For the operation of updating the model's coefficients, computations of nine points with SHM of degree 820 are needed. Because the c m , s m are independent of λ, Equation (11) needs to be evaluated only once for computation points densely spaced along a parallel (i.e., φ constant), and the efficiency can be further increased for points equally spaced in longitude [17]. According to the time complexity of a single point analyzed in our previous work, the operations of the 3 × 3 grids are approximate to (13ADD+18MUL+3DIV+1COMP)(n max + 2)(n max − 1), where n max is equal to 820. Compared with this process of gravity computation, the execution time of the fit algorithm to discrete data could be ignored. Without additional stored parameters, the space complexity of execution of coefficients update is equal to 2(n max + 4)(n max − 1) float. When n max = 820, the storage occupied is 11.4 MB. Toward a processor with 2.6 GHz master frequency, the computation time is on average approximately 2.54 s. It means that the prediction of navigation area should be finished to retain the update time for model coefficients before the vehicle reaches the boundary of the present navigation area.

Experimental Results
To illustrate the outstanding performance of PM in free-INS, a ship test was performed in the South China Sea using the Yuan Wang 5# surveying ship (shown in Figure 10). The high-precision strapdown INS used in the test consists of three ring laser gyroscopes with bias stability of 0.003 • /h and three quartz accelerometers with bias stability of 5 µg, where PM is carried out. The travel profile and DOVs on the trajectory are shown in Figures 11 and 12.

Experimental Results
To illustrate the outstanding performance of PM in free-INS, a ship test was performed in the South China Sea using the Yuan Wang 5# surveying ship (shown in Figure 10). The high-precision strapdown INS used in the test consists of three ring laser gyroscopes with bias stability of 0.003 °/h and three quartz accelerometers with bias stability of 5 µg, where PM is carried out. The travel profile and DOVs on the trajectory are shown in Figures 11 and 12.

Experimental Results
To illustrate the outstanding performance of PM in free-INS, a ship test was performed in the South China Sea using the Yuan Wang 5# surveying ship (shown in Figure 10). The high-precision strapdown INS used in the test consists of three ring laser gyroscopes with bias stability of 0.003 °/h and three quartz accelerometers with bias stability of 5 µg, where PM is carried out. The travel profile and DOVs on the trajectory are shown in Figures 11 and 12.

Experimental Results
To illustrate the outstanding performance of PM in free-INS, a ship test was performed in the South China Sea using the Yuan Wang 5# surveying ship (shown in Figure 10). The high-precision strapdown INS used in the test consists of three ring laser gyroscopes with bias stability of 0.003 °/h and three quartz accelerometers with bias stability of 5 µg, where PM is carried out. The travel profile and DOVs on the trajectory are shown in Figures 11 and 12.        According to results in Figures 13-15, the proposed polynomial model has the best performance compared with the NGM and SHM of degree 12. Relative to NGM, the maximum value of position error improvements is 500 m during the ten hours travel when navigation scheme is compensated by the polynomial model, while, it is 152 m in the navigation scheme compensated by the SHM of degree 12. Therefore, the effectiveness of the online gravity polynomial model is verified in that it could improve the position error of ship-borne INS by about ten percent.

Conclusions
The main motivation of this paper is to provide an efficient method for the application of high-degree SHM in free-INS. Due to time and space limitations, the high-degree SHM is not suitable for real-time navigation solutions. To reduce the computational complexity, we firstly propose an approximation model that converts the global support SHM to local support. Mathematical analysis illustrates the high accuracy of a low order polynomial approximating the SHM in a small area. Numerical tests in two typical areas (the North China Plain and the Himalayas)   According to results in Figures 13-15, the proposed polynomial model has the best performance compared with the NGM and SHM of degree 12. Relative to NGM, the maximum value of position error improvements is 500 m during the ten hours travel when navigation scheme is compensated by the polynomial model, while, it is 152 m in the navigation scheme compensated by the SHM of degree 12. Therefore, the effectiveness of the online gravity polynomial model is verified in that it could improve the position error of ship-borne INS by about ten percent.

Conclusions
The main motivation of this paper is to provide an efficient method for the application of high-degree SHM in free-INS. Due to time and space limitations, the high-degree SHM is not suitable for real-time navigation solutions. To reduce the computational complexity, we firstly propose an approximation model that converts the global support SHM to local support. Mathematical analysis illustrates the high accuracy of a low order polynomial approximating the SHM in a small area. Numerical tests in two typical areas (the North China Plain and the Himalayas)   According to results in Figures 13-15, the proposed polynomial model has the best performance compared with the NGM and SHM of degree 12. Relative to NGM, the maximum value of position error improvements is 500 m during the ten hours travel when navigation scheme is compensated by the polynomial model, while, it is 152 m in the navigation scheme compensated by the SHM of degree 12. Therefore, the effectiveness of the online gravity polynomial model is verified in that it could improve the position error of ship-borne INS by about ten percent.

Conclusions
The main motivation of this paper is to provide an efficient method for the application of high-degree SHM in free-INS. Due to time and space limitations, the high-degree SHM is not suitable for real-time navigation solutions. To reduce the computational complexity, we firstly propose an approximation model that converts the global support SHM to local support. Mathematical analysis illustrates the high accuracy of a low order polynomial approximating the SHM in a small area. Numerical tests in two typical areas (the North China Plain and the Himalayas) According to results in Figures 13-15, the proposed polynomial model has the best performance compared with the NGM and SHM of degree 12. Relative to NGM, the maximum value of position error improvements is 500 m during the ten hours travel when navigation scheme is compensated by the polynomial model, while, it is 152 m in the navigation scheme compensated by the SHM of degree 12. Therefore, the effectiveness of the online gravity polynomial model is verified in that it could improve the position error of ship-borne INS by about ten percent.

Conclusions
The main motivation of this paper is to provide an efficient method for the application of high-degree SHM in free-INS. Due to time and space limitations, the high-degree SHM is not suitable for real-time navigation solutions. To reduce the computational complexity, we firstly propose an approximation model that converts the global support SHM to local support. Mathematical analysis illustrates the high accuracy of a low order polynomial approximating the SHM in a small area. Numerical tests in two typical areas (the North China Plain and the Himalayas) show that the DOVs' approximation errors of a second-order polynomial are as much as 0.01 arc·s in a smooth area and 0.1 arc·s in a rugged area than full degree SHM. To accommodate the real-time navigation, we propose a novel navigation scheme assisted with an external gravity computer in which the coefficients of the polynomial model will update synchronously. Further studies demonstrate that the optimal degree of basic SHM and the number of grid points for fitted coefficients are 820 and 3 × 3, respectively. Combined with the theoretical analysis, the proposed gravity model is finally expressed as a second-order polynomial in each 6 × 6 navigation area. Computational complexity and largest l2-norm error of the model are very low as compared with the traditional NGM, typical modified NGM, the cubed-sphere model, and SHM of degree 12 and 820. The efficiency of polynomial model is verified in a ship test. Compared with traditional NGM, the maximum position errors in the travel profile were decreased by ten percent.