The content of this article is from the 8th issue of the Journal of Surveying and Mapping in 2024 (drawing review number: GS Jing (2024) No. 1748)
A ico_HEALPix grid-based modeling method for ultra-high-order earth gravity field
Zhang Zhanpeng
,1, Li Xinxing,1,2,3, Liu Changjian1, Fan Haopeng1, Pei Xianyong1
1. School of Geospatial Information, Information Engineering University, Zhengzhou 450001, Henan, China
2. State Key Laboratory of Geographic Information Engineering, Xi'an 710054, Shaanxi, China
3. Key Laboratory of Smart Earth, Beijing 100020, China
Summary:
National Natural Science Foundation of China (U23A2028) (41404020) (42174007) (42074013); State Key Laboratory of Geographic Information Engineering(SKLGIE2023-Z-1-1)
About the Author:
张展鹏（2000—），男，硕士生，研究方向为测量数据处理理论与方法。 E-mail:[email protected] 通讯作者: 李新星 E-mail:[email protected]
Summary:
In order to solve the problem of data redundancy in the construction of the gravity field model segmented by traditional geographic grid data in the high latitude region, this paper introduces the hierarchical equal area equal latitude pixelated (HEALPix) grid structure into the calculation of the earth's gravity field, and proposes an ultra-high-order earth gravity field modeling method using icosahedral HEALPix(ico_HEALPix) grid, which realizes the efficient construction of the global 3600-order subspherical harmonic coefficient, and at the same time, it is aimed at ico_ In the process of spherical harmonic analysis, the method matrix of HEALPix mesh is not a problem of strict block diagonalization structure, and an iterative algorithm is designed to effectively improve the accuracy of model construction. The experimental results show that under the premise that the data volume of the ico_HEALPix grid data is less than 5 million in the geographic grid, the accuracy of the global earth gravity field model constructed by the iterative method can be better than that of the geographic grid, and the error order RMS of the spherical harmonic coefficient is increased by 1~2 orders of magnitude, and the problems of distortion and data redundancy at the north and south poles of the geographic grid are solved, and the data utilization rate of the grid is improved.
keyword
HEALPix grids; ico_HEALPix grid; ultra-high-order Earth gravity field model; iterative algorithms; least squares; XGM2019e
This article cites format
Zhang Zhanpeng, Li Xinxing, Liu Changjian, Fan Haopeng, Pei Xianyong. Ultra-high-order Earth Gravity Field Modeling Method Based on ico_HEALPix Mesh[J]. Journal of Surveying and Mapping, 2024, 53(8): 1531-1539 doi:10.11947/j.AGCS.2024.20230151
ZHANG Zhanpeng, LI Xinxing, LIU Changjian, FAN Haopeng, PEI Xianyong. Ultra-high-order Earth gravity field modeling method based on ico_HEALPix grid[J]. Acta Geodaetica et Cartographica Sinica, 2024, 53(8): 1531-1539 doi:10.11947/j.AGCS.2024.20230151
Read more
http://xb.chinasmp.com/article/2024/1001-1595/1001-1595-2024-08-1531.shtml
The Earth gravity field model based on the traditional conformal geographic grid faces many insurmountable limitations, including unequal grid area, redundant observation data, complex gravity field smoothing factor, and large integral discretization error. In view of the problems faced by the traditional geographic grid in spherical harmonic analysis and solving, many scholars have adopted a new meshing method to solve the relevant problems and achieved good results. Ref. [6] shows that the quality of the normal matrix structure for solving the Earth's gravity field model is affected by the spherical mesh distribution, and the spherical uniformly distributed mesh structure is very important for spherical harmonic analysis and synthesis problems. In Ref. [7], the fast fourier transform (FFT) technique was used to verify that the approximate global uniform distribution of points based on octahedra and icosahedral has significantly improved efficiency and accuracy in spherical harmonic analysis and synthesis compared with the geographic grid. In astronomy, analyzing the cosmic microwave background (CMB) data of the celestial sphere can better understand the geometry and composition of the universe [8], and in order to solve the needs of visualization and rapid data analysis of discrete data on the celestial sphere, hierarchical equal area isolatitude pixelation (SIPIXELATION) has been designed in Ref. [9-10]. HEALPix) grid structure, due to the excellent characteristics of the grid and the release of the software of the same name [11], the HEALPix partitioning scheme has been popularized in applications such as astrophysics [12], astronomy [13-14], geophysics [15], planetary science [16], nuclear engineering [17], and computer vision [18]. In Ref. [19], a new combination of inhomogeneous FFT, double Fourier spherical method, and Slevinsky fast spherical harmonic transform was used to perform spherical harmonic analysis of a grid type based on HEALPix, which showed that the convergence speed was at least two times higher while improving the accuracy of spherical harmonic analysis. In order to solve the problems of uneven data distribution and redundancy of high-latitude data in the construction of the earth's gravity field model from the gravity data of the geographic grid, this paper introduces the HEALPix meshing scheme into the spherical harmonic analysis, which proves that the conventional HEALPix mesh has a large error due to insufficient sampling in the longitude direction of the high-latitude data in the spherical harmonic analysis ico_HEALPix. least-squares, LS) [20-22], an iterative algorithm is designed to solve the problem of incomplete block diagonalization of normal matrices, and the construction of a 3600-order earth gravity field model is realized.
1 Modeling Methodology
1.1 HEALPix Grid and its Promotion
The HEALPix grid system is a multi-level composite grid model that uses spherical triangular projections in polar regions and cylindrical equal-area projections in non-polar regions. It has three core features [23]:(1) the data is constructed hierarchically according to resolution; (2) the mesh area of the spherical surface is equal; (3) The grid is distributed at equal latitude on the spherical surface. The spatial and planar grids of the model are shown in Figure 1.
Figure 1
Fig.1 Spatial grids [19] and planar grids
Fig.1 Spatial grid[19]and planar grid
The HEALPix mesh system consists of 12 basic grids, divided into 3 rows from north to south, each row is evenly distributed with 4 basic grids, each basic grid can subdivide the edge length into Nside=2k parts, where k=1,2,3,...,n represents the level of meshing, so each basic grid is divided into
Subgrids of equal area of blocks. The ico_HEALPix (icosahedral HEALPix) mesh is a modification of the conventional HEALPix mesh in this paper, which expands the 12 basic grids of HEALPix mesh 3×4 to a total of 20 basic grids × 45, and its plane expansion is shown in Figure 2(b). There are 4 rows from north to south, and each row is evenly distributed with 5 basic grids. In order to facilitate the calculation of the gravity field, the overall translation of the southern hemisphere of the 4×5 basic grid is π/5 to obtain a symmetrical grid about the north-south of the equator, and the plane distribution of the translated grid is shown in Figure 2(c). ico_HEALPix retains the excellent characteristics of the HEALPix grid system, such as hierarchical construction according to resolution and equal latitude distribution of grid area.
Figure 2
Fig.2 HEALPix mesh and its improved mesh plane deployment
Fig.2 HEALPix grid and its improved grid plane expansion
As shown in Figure 3, the planar expansion and global distribution map of the ico_HEALPix grid corresponding to k=1 are indicated. The edges represent the boundaries of the mesh, and the solid points represent the location of the center point for each small geographic grid.
Figure 3
Fig. 3 ico_HEALPix grid plane unfolding and global distribution
Fig.3 Plane expansion and global distribution of ico_HEALPix grid
Eq. (1) gives the longitude resolution λi corresponding to different latitude circles i from north to south
(1)
ico_HEALPix grids have subgrids of equal size but different shapes
blocks, the center points of these grids are located on 5×Nside-1 latitude circles, and the pixel size of each grid per unit sphere is
, each pixel can calculate its coordinates (θ, λ) on the spherical surface using the index value of the grid latitude circle (Nside, i, j), where i represents the ith layer of latitude circle from north to south, j represents the column number where the corresponding latitude circle is located, θ is the coincidence of the sphere, 0 at the North Pole, π at the South Pole, λ is the longitude, and the range [0, 2π]. The Arctic region is defined when 1≤i < Nside, the Antarctic region when 4Nside <i<5Nside, and the equatorial region between the Arctic region and the Antarctic region, at which point Nside ≤ i≤4Nside. Let z=cos θ, and the formula for calculating the coordinates of the center point of the grid is as follows. In the Arctic region, the index value (i,j) satisfies 1≤i<Nside, 1≤j≤5i
(2)
(3)
(4)
where s denotes the offset of the index value j relative to latitude 0 in the latitude circle index system.
赤道区域，索引值(i,j)满足Nside≤i≤4Nside,1≤j≤5Nside
(5)
(6)
(7)
南极区域，索引值(i,j)满足4Nside<i<5Nside,1≤j≤25Nside-5i
(8)
(9)
(10)
With the above formula, the spherical coordinates (θ, λ) of the center point of the grid on the sphere can be calculated from the current index value (i,j), and the value of s is 0 or 1.
In order to fully illustrate the advantages of the ico_HEALPix, the structural shapes of the two grids are compared under the condition that the number of geographic grids and ico_HEALPix grids is roughly equal. The geographic grid with a resolution of 18°×18° and the ico_HEALPix grid with k=1 were selected, and the comparison of the number of points between the two grids is shown in Table 1.
Table 1 Comparison of the number of geographic grids and ico_HEALPix grids
Tab.1 Comparison of the number of geographic grids and ico_HEALPix grids
| grid | Use model orders | Number of meshes |
|---|---|---|
| Geographic grids | 7 | 98 |
| ico_HEALPix grid | 7 | 80 |
A new window opens| Download the CSV
Figure 4 shows the distribution of the Earth's gravitational field in the Arctic with the geogrid and the improved HEALPix grid.
Figure 4
Fig.4 Comparison of the structure of the Arctic region between the geographic grid and the ico_HEALPix grid
Fig.4 Comparison of Arctic structure between geographic gird and ico_HEALPix grid
As shown in Figure 4(a), the area of the geographic grid based on latitude and longitude division is not equal from the poles to the equator, and the grid area of different latitude circles varies greatly, or even not in the same magnitude, and there is a large distortion problem. As shown in Figure 4(b), the area distortion problem in the Arctic and Antarctic regions has been effectively solved ico_HEALPix the equal latitude distribution of the grid, and the grid is evenly and seamlessly distributed. In order to illustrate the advantages of ico_HEALPix grid over geographic grids in the equatorial region, as shown in Figure 5, the grid ranges from 17°E to 110°E and 30°N to 30°S were intercepted for comparative analysis.
Figure 5
Fig.5 Comparison of equatorial regions between geographic grids and ico_HEALPix grids
Fig.5 Comparison of geographic grid and ico_HEALPix grid equatorial region
As can be seen from Figure 5, the shape of the equatorial region of the geographic grid is too large, and the ico_HEALPix grid of equal area partition makes the grid area of the equatorial region smaller than that of the geographic grid, and the resolution of the equatorial region is improved. The overall structure shape is shown in Figure 6.
Figure 6
Fig.6 Comparison of the overall structure of the geographic grid and the ico_HEALPix grid
Fig.6 Comparison of the overall structure between geographic grid and ico_HEALPix gird
Under the condition that the number of grids is approximately equal, compared with the geographic grid, the ico_HEALPix grid has the following three advantages: (1) the equal area segmentation of the ico_HEALPix grid solves the problem of over-density and low data utilization rate of the geographic grid; (2) The equal area partitioning of the ico_HEALPix grid makes the data resolution of the grid better than that of the geographic grid in the equatorial region, and the grid structure is more reasonable. (3) The ico_HEALPix mesh is approximately evenly distributed on the spherical surface, and the isotropy is stronger than that of the geographic grid.
1.2 The relationship between gravity anomalies and potential coefficients
The spherical harmonic expansion formula for gravity anomalies is [24-25]
(11)
where Δg is the observation of the spatial gravity anomaly; (r,θ,λ) denotes the spherical coordinates of the gravitational anomaly;
is the perturbation potential coefficient of the gravity field model for the obtained set of gravity field models, i.e., the quantity to be found; a is the long radius of a normal ellipsoid;
(cos θ) is the complete normalized association Legendre function; DEG is the highest order of the spherical harmonic model. Eq. (11) is the gravitational outlier Δg and the bit coefficient
The bridge between the two is also the basis for solving the bit coefficient by least squares. Expresses observational equation (11) as a least-squares generic expression
(12)
where,
(13)
The F-function is a linear equation with respect to the spherical harmonic coefficient Xa, and according to the least-squares adjustment model, there is a legal equation [26]
(14)
where N=ATPA is called the coefficient matrix of the normal equation, which is referred to as the normal matrix; U=ATPLb is the vector of the free term of the normal equation. The individual elements of matrix A are calculated as:
(15)
(16)
where (θi,λj) is the colatitude and longitude of the center of the sphere corresponding to the index value (i,j) of the gravity anomaly grid data. Its solution is:
(17)
The above solution process needs to satisfy the number of observations of the gravity anomaly Δg is greater than the spherical harmonic coefficient to be calculated
For ico_HEALPix and HEALPix meshes, only the total number of grid points needs to be satisfied with the number of coefficients to be found, while for the geographic grid distribution, dθ=dλ needs to satisfy the Nyquist sampling law, i.e., Deg≤π/dθ, and the equal sign cannot be taken when using the least squares method, i.e., Deg<π/dθ [21]. Table 2 shows the recoverable model orders of ico_HEALPix mesh at different resolutions and the minimum number of geogrid points required for the corresponding order model.
Table 2 Comparison of the number of points required between grids
Tab.2 Comparison of the required number of points between grids
| You | k | Nside | The number of bit coefficients | ico_HEALPix points | HEALPix Credits | Number of geogrid points |
|---|---|---|---|---|---|---|
| 180 | 6 | 64 | 32 757 | 81 920 | 49 152 | 64 800 |
| 360 | 7 | 128 | 130 317 | 327 680 | 196 608 | 259 200 |
| 720 | 8 | 256 | 519 837 | 1 310 720 | 786 432 | 1 036 800 |
| 1440 | 9 | 512 | 2 076 477 | 5 242 880 | 3 145 728 | 4 147 200 |
| 1800 | 9 | 512 | 3 243 597 | 5 242 880 | — | 6 480 000 |
| 1800 | 10 | 1024 | 3 243 597 | — | 12 582 912 | 6 480 000 |
| 2160 | 10 | 1024 | 4 669 917 | 20 971 520 | 12 582 912 | 9 331 200 |
| 3600 | 10 | 1024 | 12 967 197 | 20 971 520 | — | 25 920 000 |
| 3600 | 11 | 2048 | 12 967 197 | — | 50 331 648 | 25 920 000 |
| 5400 | 11 | 2048 | 29 170 797 | 83 886 080 | 50 331 648 | 58 320 000 |
Note: The "—" in Table 2 indicates that the number of grid points is less than the number of bit coefficients, which cannot be used to solve the spherical harmonic coefficients of the corresponding order model.
A new window opens| Download the CSV
In order to illustrate the advantages of the proposed ico_HEALPix mesh in calculating the spherical harmonic coefficient, the HEALPix mesh data of k=11 and the ico_HEALPix mesh data of k=10 were selected to recover the bit coefficients of the 5400th and 3600th order models, respectively.
1.3 Block diagonal least-squares iteration method
XGM2019e is a set of 5540-order combined global gravity field models published by the Technical University of Munich in 2019, corresponding to a spatial resolution of 2' [27]. In this paper, the first 6th order model coefficients of XGM2019e are used to calculate the gravity outliers of the points in the ico_HEALPix grid with k=1, and then the least squares method is used to solve the spherical harmonic coefficients, and the N coefficients of the normal matrix are shown in Figure 7, where the white part represents the coefficient value close to zero, and the black part represents the non-zero value, which shows that the normal matrix of the 6th order earth gravity field model based on the ico_HEALPix grid distribution is a non-complete block diagonal matrix, and there are non-zero values outside the block diagonal elements. Therefore, it is difficult to solve the large normal inverse matrix N-1.
Figure 7
Fig.7 ico_HEALPix grid matrices of k=1
Fig.7 The ico_HEALPix grid method matrix of k=1
Because N is a symmetric positive definite matrix with dominant diagonals, the "blocks" in the normal matrix N are inverted one by one, ignoring the non-block diagonal matrix parameters, and an approximation N* to N-1 can be obtained, and the cyclic iterative solution formula can be obtained by using this matrix
(18)
where r(k+1) is the residual vector of the k+1st iteration. Combine Eq. (18) to get a new expression
(19)
When the spectral radius of the matrix (I-N*N) is less than 1, the above iterative method converges, and the final iteration result is consistently approximate to the true value.
To solve this problem, an iterative algorithm is designed, and the specific iteration process is shown in Figure 8.
FIGURE 8
Figure 8 Diagonal least squares iteration process of blocks
Fig.8 Block diagonal least squares iterative flow chart
2. Numerical simulation test and analysis
2.1 Comparison of HEALPix grids with traditional geographic grids
In this paper, we propose to ico_HEALPix the mesh instead of using the HEALPix mesh because of the limitations of using the HEALPix distribution to sample gravity anomalies and recover the Earth's gravity field model. In order to illustrate this point, the HEALPix grid distribution of k=6 was used as the sampling point, the gravity outliers of the model simulated by the model 180-order coefficients were XGM2019 at the sampling points, the 180-order model bit coefficients were recovered by the block diagonal least squares method, and then the gravity anomalies of the model points at the sampling points were positively calculated, and the difference distribution was shown in Figure 9 compared with the original simulation values.
Figure 9
Fig.9 HEALPix mesh recovers the difference caused by the spherical harmonic coefficient
Fig.9 The difference caused by HEALPix grid restoring spherical harmonic coefficient
As can be seen from Figure 9, the gravity anomaly of the proposed model is in good agreement with the original simulation value in the middle and low latitude regions. However, in the high latitude region, there are several mGal differences between the gravity anomaly of the proposed model and the original simulation value. Further analysis suggests that this difference is due to the fact that the sampling of HEALPix mesh in the longitude direction of the high-latitude region is too sparse. In order to fully illustrate the problems existing in HEALPix mesh, the model coefficients of the first 5400 orders of XGM2019e were selected as the input data, and the gravity outliers of the center points of the geographic grid with a resolution of 2'×2' and the HEALPix grid with a partition level of k=11 were simulated, and the required number of grid points is shown in Table 3.
Table 3 Comparison of the number of geographic grids and HEALPix grids
Tab.3 Comparison of the number of geographic grids and HEALPix grids
| grid | Order | Number of meshes |
|---|---|---|
| Geographic grids | 5400 | 58 320 000 |
| HEALPix grids | 5400 | 50 331 648 |
A new window opens| Download the CSV
The gravity anomaly data of the center point of the geographic grid and the HEALPix grid were obtained by simulation, and then the model coefficients were backcalculated by using the block diagonal least squares iterative method, and compared with the true values of the XGM2019e model, the order error RMS of the calculated model bit coefficients was calculated, and the calculation formula was [28]
(20)
The process shown in Figure 8 is used to iteratively calculate the ultra-high-order sub-earth gravity field model constructed by HEALPix grid, and the order error RMS of non-iteration, 4 iterations, and 9 iterations are counted, and the accuracy results of the model solved by the geographic grid are compared, as shown in Figure 10. Because the geogrid matrices strictly satisfy the diagonal form of blocks, there is no iteration.
FIGURE 10
Fig.10 Comparison of error-order RMS
Fig.10 Error order RMS comparison
The error spectra of the HEALPix grid without iteration, the HEALPix grid iteration 9 times, and the geographic grid are shown in Figure 11.
FIGURE 11
Fig.11. Error spectrum of spherical harmonic coefficient of the model
Fig.11 Model spherical harmonic coefficient error spectrum
As can be seen from Fig. 10 and Fig. 11, the low-order accuracy of the earth gravity field model based on HEALPix mesh has been improved to a certain extent, but the accuracy of the high-order part of the model has been poorly restored, which is due to the low sampling rate of the HEALPix mesh in the longitude direction of the high-latitude region.
2.2 Comparison of ico_HEALPix with traditional geographic grids
The sampling rate of the ico_HEALPix grid in the longitude direction is 1.25 times that of HEALPix, which can better improve the problem that the HEALPix grid is too low in the longitude direction in the high latitude region, and the simulation experiment is used to verify it. The model coefficients of the first 3600 orders of XGM2019e with ico_HEALPix grid points less than the geographic grid were selected as the input data, and the gravity outliers at the points in the geographic grid with a resolution of 3'×3' and the ico_HEALPix grid with a partitioning level k=10 were simulated as the observation data, and the model bit coefficients of the first 3600 orders were recovered by using the block diagonal least squares iteration method. Using XGM2019e as the standard model, the geo-grid and ico_HEALPix-grid recovered the model bit coefficients of 3600 orders, respectively, and the required number of meshes is shown in Table 4.
Table 4 Comparison of the number of ico_HEALPix grids and geographic grids
Tab.4 Comparison of the number of geographic grids and ico_HEALPix grids
| grid | Order | Number of meshes |
|---|---|---|
| Geographic grids | 3600 | 25 920 000 |
| ico_HEALPix grid | 3600 | 20 971 520 |
A new window opens| Download the CSV
The process shown in Figure 8 is used to iteratively calculate the ultra-high-order sub-earth gravity field model constructed ico_HEALPix the grid, and the order error RMS of the non-iterative approximate solution, 4 iterations, and 9 iterations are counted, and the model accuracy results are compared with the model accuracy results of the geographic grid solution, as shown in Figure 12.
FIGURE 12
Fig.12 Comparison of error order RMS
Fig.12 Error order RMS comparison
The error spectra for ico_HEALPix grid no iteration, ico_HEALPix grid iteration 9 times, and geographic grid are shown in Figure 13.
FIGURE 13
Fig.13. Error spectrum of spherical harmonic coefficient of the model
Fig.13 Model spherical harmonic coefficient error spectrum
As can be seen from Fig. 12 and Fig. 13, the accuracy of the global earth gravity field model constructed by the iterative method can be better than that of the geographic grid under the premise that the data volume of the ico_HEALPix grid data is less than that of the geographic grid by about 5 million, and the coefficient error order RMS is increased by 1~2 orders of magnitude.
3 Conclusion
In view of the limitations of geographic meshing in the process of spherical harmonic analysis of gravity field, including the obvious mixing effect in the process of model solving, the limitation of solution accuracy, and the inherent contradiction between the Nyquist sampling criterion and the data resolution, this paper introduces spherical equal-area ico_HEALPix meshing into the model solving, and designs a block least squares iterative algorithm, which realizes the high-precision construction of an ultra-high-order earth gravity field model based on ico_HEALPix mesh at the cost of losing a certain computational efficiency. Taking the XGM2019e model as the input truth, the results of spherical harmonic synthesis and spherical harmonic analysis show that the accuracy of the ultra-high-order subterrestrial earth gravity field model based on the ico_HEALPix grid can be better than that of the geographic grid after iteration, and at the same time, the scheme improves the resolution accuracy of the data that can be stored in the original model, reduces the data storage capacity under the same circumstances, solves the distortion problem near the north and south poles, and improves the utilization and reliability of the grid data. In the future, other new uniform meshing methods can be explored to avoid the iterative solution of the normal matrix, so as to improve the calculation efficiency and ico_HEALPix the application advantages of the mesh in weakening the mixing effect. The model constructed in this paper can provide reference for improving the accuracy of spherical harmonic coefficient model, enriching the theory of constructing the earth's gravity field model, and improving the theoretical system of spherical harmonic coefficient model.
First trial: Zhang Yanling review: Song Qifan
Final Judge: Jin Jun
information
○ Professor Shi Yan, Central South University: Spatiotemporal Anomaly Detection: Connotation Transformation from Data-Driven to Knowledge-driven and Realization Path|Journal of Surveying and Mapping, Vol. 53, No. 8, 2024