A grid database generation and indexing method for gravity compensation interpolation calculation

Abstract

The application provides a grid database generation and index method for gravity compensation interpolation calculation.The method comprises the following steps: generating a global gravity grid by using a high-precision global gravity model; considering the spherical change of the earth, appropriately cutting the number of gravity grid points to form a gravity grid with uniform spherical distribution; determining the data minimum significant digit according to the demand of an inertial reference system, and storing the uniform gravity grid points in binary form to form a uniform gravity grid database; and indexing the required gravity grid data through the pre-stored corner point information and the position information of the interpolation calculation point.The application can be applied to a gravity compensation module of a high-precision inertial reference system, and the method can not only reduce the database storage space, but also ensure the accuracy of the gravity compensation result and the index efficiency.Meanwhile, the method can flexibly adjust the generation mode of the gravity grid database according to the demand of the inertial reference system, and has the characteristics of simple realization and strong engineering applicability.

Description

Technical Field

[0001] This invention belongs to the field of inertial navigation technology, and specifically relates to a grid database generation and indexing method for gravity compensation interpolation calculation. Background Technology

[0002] In inertial navigation calculations, the accuracy of gravity compensation directly affects the accuracy of the output information of the inertial reference system. Currently, gravity compensation in inertial reference system calculation algorithms mostly adopts the normal gravity field model, which describes the Earth's gravity field under the assumption of a regular Earth. However, in reality, the Earth's mass distribution is non-uniform, and the normal gravity field can only approximately represent the Earth's true gravity field. Accurately calculating the Earth's true gravity field is of great significance for high-precision inertial navigation. Compared to the normal gravity field model, a high-precision gravity field model can improve inertial navigation accuracy by approximately 0.1 n mile / h; simultaneously, a 1 mGal horizontal gravity disturbance can cause approximately 0.2 arcseconds of horizontal attitude error (horizontal gravity disturbances are typically tens to 100 mGal).

[0003] For high-precision inertial reference systems, a gravity compensation module is needed to accurately calculate the Earth's true gravity field in order to improve the system's output accuracy. Currently, the most accurate gravity compensation scheme is to directly calculate the gravity value at a point using higher-order spherical harmonics, but this method requires extremely high computational resources. Considering factors such as cost, existing civilian inertial reference systems cannot meet these computational requirements. Therefore, considering engineering applications, interpolation calculation schemes have become the preferred solution in the engineering field.

[0004] Although interpolation calculation schemes can reduce the computational burden of the computing system, the databases used by the interpolation calculation schemes will increase the storage burden of the system. To reduce the computational burden of the system, a grid database with low storage capacity needs to be designed. Summary of the Invention

[0005] The purpose of this invention is to address the aforementioned technical problems in the prior art by providing a grid database generation and indexing method for gravity compensation interpolation calculation. In order to reduce the storage burden of the system, this method can generate a gravity grid with uniform spherical distribution.

[0006] The technical solution of this invention: This invention provides a grid database generation and indexing method for gravity compensation interpolation calculation. Global gravity grid data is generated using a high-precision global gravity model; based on the Earth's spherical surface variations, the generated global gravity grid is trimmed to produce gravity grid data with a uniform spherical distribution; the least significant bit of the data is determined according to the requirements of the inertial reference system, and the uniformly distributed gravity grid is stored in binary form to form a uniform gravity grid database; the required gravity grid data is indexed using pre-stored corner information and the position information of interpolation calculation points.

[0007] A global gravity grid is a latitude-longitude grid that stores gravity data at fixed longitude and latitude intervals worldwide. Global gravity grid data represents the gravity data at the intersections of meridians and parallels on the global gravity grid.

[0008] A uniformly distributed spherical gravity grid is a latitude-longitude grid that stores gravity data at different latitudinal intervals across different latitudinal zones. The spherical distances between the intersections of meridians and parallels are not significantly different from the spherical distances between the intersections of meridians and parallels at the equator. Uniformly distributed gravity grid data represents gravity data at the intersections of meridians and parallels on a globally distributed spherically distributed gravity grid.

[0009] Specifically, the process includes the following steps: S1: Generate global gravity grid data using a high-precision global gravity model; S2: Based on the Earth's spherical variations, trim the latitudinal directions of the high-latitude grid in the generated global gravity grid data to generate gravity grid data with a uniform spherical distribution; S3: Determine the least significant bit of the data according to the requirements of the inertial reference system, and store the uniformly distributed gravity grid data in binary form to form a uniform gravity grid database; S4: Divide the data into latitudinal data blocks and store the information of each data block; S5: Index the data block where the carrier is located based on the information of each data block and the carrier's position coordinates; S6: Calculate the interpolation grid point data required to determine the gravity value at the carrier's location on the data block where the carrier is located.

[0010] In one possible embodiment, the process of generating global gravity grid data in step S1 specifically includes: selecting a latitude and longitude grid with a fixed interval of δ arcminutes, and using the global ultra-high-order Earth gravity field model EGM-2008 to calculate and record the gravity field data at the intersection of latitude and longitude lines.

[0011] In one possible embodiment, in step S2, for regions with latitude |θ|>60°, the global gravity grid is cropped, and the longitude interval of the gravity database is adjusted. Adjusted to:

[0012]

[0013] Here, θ refers to latitude.

[0014] In one possible embodiment, in step 3, the least significant bit of data is determined based on the requirements of the inertial reference system, where the requirements of the inertial reference system refer to the zero-bias instability of the accelerometer, and the least significant bit refers to the number of significant digits that need to be retained.

[0015] High-precision gravity values ​​([GE, GN, GU]) are calculated using high-precision spherical harmonic functions. The celestial component GU of the high-precision gravity value is then subtracted from the celestial component GU_Normal calculated by the normal gravity formula to obtain the celestial residual component.

[0016] Determine the zero-bias instability Kμg (in μg) of the accelerometer in the inertial reference system;

[0017] One-tenth of the accelerometer's zero-bias instability is used as the basis for determining the effective bit, K / 10μg;

[0018] Since high-precision gravity values ​​calculated using high-precision spherical harmonic functions are often expressed in grams, and K / 10μg = K × 10^(-7)g, grid points need to be retained to at least Q decimal places (Q = (7 - num(K) + 1)). (It is stipulated that K < 1000. If K is greater than 1 and less than 10, then num(K) = 1; if K is greater than or equal to 10 and less than 100, then num(K) = 2; if K is greater than or equal to 100 and less than 1000, then num(K) = 3.)

[0019] Multiply the [GE, GN, GU] components by 10^Q to obtain integers, discard the decimal parts, and then convert the integers into binary for storage.

[0020] In one possible embodiment, in step S4, for latitude zones with latitude |θ|>60°, they are symmetrically divided into 2N_1 data blocks along the equatorial plane; together with the data blocks for latitude zones with latitude |θ|<60°, a total of 2N_1+1 data blocks are formed; the data block information includes the upper latitude boundary, the lower latitude boundary, the number of grid points in the data block, and the longitude interval.

[0021] In one possible embodiment, in step S5, the coordinates of the carrier's position are [λ φ h], where λ is latitude in degrees; φ is longitude in degrees; and h is altitude in meters.

[0022] The upper and lower latitude boundaries of each data block information obtained in step S4 are denoted as the dataset CORNER_BLOCK[i]. k上 i k下 Where i k上 Indicates the upper dimensional boundary of the k-th data block; i k下 The upper dimension represents the lower boundary of the k-th data block; k = 1, ..., 2N_1+1; i k上 =i (k-1)下 i k下 =i (k+1)上 i 1上 =90; i (2N_1+1)下=-90; The number of grid points in each data block is denoted as dataset DATA_num=[NN k ], k = 1, ..., 2N_1+1; where N ∈ N k This represents the number of grid points in the k-th data block;

[0023] Let k = N_1 + 1;

[0024] If λ-i k下 If >0, then continue calculating λ-i. k上 If λ-i k上 If λ-i < 0, then the carrier is determined to be on this data block; if λ-i k上 If the value is greater than 0, then continue setting k = k - 1 until λ - i is found. k上 <0;

[0025] If λ-i k下 If λ < 0, then let k = k + 1 and continue calculating λ - i. k下 If λ-i k下 If >0, then calculate λ-i k上 If λ-i k上 If λ-i < 0, then the carrier is determined to be on this database block; if λ-i k上 If the value is greater than 0, continue setting k = k + 1 until λ - i is found. k上 <0; Determine the data block number k where the carrier is located, denoted as k1.

[0026] In one possible embodiment, in step S6, let the larger absolute value of the upper and lower latitude boundaries of the data block where the carrier is located be CORNER_BLOCK_MAX, and the smaller absolute value be CORNER_BLOCK_MIN, and calculate |CORNER_BLOCK_MAX-λ|, where || represents the absolute value of the data;

[0027] Calculate the latitudinal interval Δλ, Δλ = floor(|CORNER_BLOCK_MAX-λ|×60 / δ), where floor represents rounding down;

[0028] Calculate the longitude interval Δφ, Δφ = floor(|φ - LAT) REQ |×60×cos(CORNER_BLOCK_MIN) / δ);

[0029] Among them, LAT REQ The starting longitude of the data block where the carrier is located;

[0030] The interpolation grid points are calculated as follows:

[0031] po int 1=NN1+……+NN k1-1 +Δφ×N+Δλ

[0032] po int 2=NN1+……+NN k1-1 +Δφ×N+Δλ+1

[0033] po int 3=NN1+……+NN k1-1 +(Δφ+1)×N+Δλ+

[0034] po int 4=NN1+……+NN k1-1 +(Δφ+1)×N+Δλ+1

[0035] Where N is (|CORNER_BLOCK_MAX-CORNER_BLOCK_MIN|)×60 / δ.

[0036] According to a second aspect of the present invention, a computer-readable storage medium is provided, the computer-readable storage medium storing a computer program, characterized in that the computer program, when executed by a processor, implements the steps of the above-described method.

[0037] The advantages of this invention are:

[0038] 1. It can reduce the storage burden on the system;

[0039] 2. It can ensure the high efficiency of gravity compensation;

[0040] 3. It can ensure that the database information required for interpolation calculations can be found quickly;

[0041] 4. Low cost and high versatility: This solution can be used for the design of gravity compensation modules in high-precision inertial reference systems. Attached Figure Description

[0042] To more clearly illustrate the technical solutions implemented in this invention, a simple explanation of the accompanying drawings used in the description of this invention will be provided below. Obviously, the drawings described below are merely some embodiments of this invention. For those skilled in the art, other drawings can be obtained based on these drawings without any creative effort.

[0043] Figure 1 A schematic diagram of the global gravity grid;

[0044] Figure 2 This is a schematic diagram illustrating the principle of gravity grid clipping.

[0045] Figure 3 A schematic diagram of a gravity grid with uniform distribution on a spherical surface;

[0046] Figure 4 This is a schematic diagram showing the carrier's position, corner points, and interpolation points;

[0047] Figure 5 A schematic diagram of the ideal interpolation calculation mesh (the interpolation algorithm assumes that the gravity anomaly and the vertical deviation are within...). Figure 5 The grid lines shown are distributed in a regular pattern.

[0048] Figure 6 This is a schematic diagram of the global gravity grid in high-latitude regions (when a uniform latitude and longitude grid is selected for interpolation calculation, the longitude grid lines will be distorted in high-latitude regions, and the interpolation effect will decrease).

[0049] Figure 7 This is a schematic diagram of the gravity grid formed in Example 1;

[0050] Figure 8 This is a schematic diagram of planar interpolation;

[0051] Figure 9 This is a schematic diagram of three-dimensional interpolation. Detailed Implementation

[0052] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0053] The features and illustrative embodiments of various aspects of the present invention will now be described in detail. In the following detailed description, numerous specific details are set forth in order to provide a thorough understanding of the invention. However, it will be apparent to those skilled in the art that the invention may be practiced without requiring some of these specific details. The following description of embodiments is merely intended to provide a better understanding of the invention by illustrating examples of the invention. The invention is by no means limited to any specific setup and method set forth below, but covers any improvements, substitutions, and modifications to the structures, methods, and devices without departing from the spirit of the invention.

[0054] It should be noted that, unless otherwise specified, the embodiments of the present invention and the features thereof can be combined with each other, and the various embodiments can be referenced and cited in each other. The present invention will now be described in detail with reference to the accompanying drawings and embodiments.

[0055] A method for generating and indexing a grid database for gravity compensation interpolation calculations includes the following steps:

[0056] Step 101: Generate global gravity grid data using a high-precision global gravity model;

[0057] A latitude and longitude grid with a longitude and latitude interval δ of 5 arcminutes was selected. Using EGM-2008, the gravity field data at the intersection points of the latitude and longitude lines were calculated and recorded. The generated gravity grid is shown below. Figure 1 As shown, the black dots represent gravity data from the global gravity grid.

[0058] Step 102: Based on the changes in the Earth's spherical surface, the generated global gravity grid is trimmed to generate gravity grid data with a uniform spherical distribution;

[0059] Considering the variations in the Earth's spherical surface, there is a significant difference between the spherical distance of 5 arcminutes near the equator and the spherical distance of 5 arcminutes near the poles. Therefore, using a 5 arcminute latitude and longitude grid for the grid near the poles would increase the storage burden considerably. Thus, it is necessary to trim the global gravity grid. The basis for this trimming will be explained in detail below.

[0060] Suppose there are two gravity grid data points on the equatorial plane, separated by an arcminute distance δ. Their corresponding spherical distance Δr can be expressed as:

[0061] Δr=δ·r e

[0062] Where r e The radius is the Earth's radius.

[0063] On the plane at latitude θ°, such as Figure 2 As shown, to obtain the longitude difference corresponding to the spherical distance Δr... The following calculation formula needs to be used:

[0064]

[0065] According to the above formula, on a plane at latitude θ° north, two points do not need to use the same longitude difference as the equatorial plane to ensure equal spherical distance. If we assume that the global gravity grid is clipped in regions where latitude |θ|>60°, and if N_1 is 1, it will form 3 data blocks with latitude zones of 90°, 60°, -60°, and -90°. Then, in regions where latitude |θ|>60°, the interval of the gravity database... It can be represented as:

[0066]

[0067] The generated spherically uniform gravity grid is as follows: Figure 3 As shown.

[0068] Step 103: Determine the least significant bit of the data according to the requirements of the inertial reference system, and store the uniformly distributed gravity grid on the spherical surface in binary form to form a uniform gravity grid database;

[0069] If the zero-bias repeatability of the accelerometer in the inertial reference system is 20 μg (1σ), to ensure the longitude of gravity compensation, the gravity grid data with uniform spherical distribution must be accurate to seven decimal places. If the gravity value at a certain point is 9.81563546655, it should be stored as 9.8156355. Furthermore, to save storage space, this gravity value also needs to be stored in binary form, thus forming a gravity grid database with uniform spherical distribution.

[0070] This invention provides an indexing method for grid points in a spherically distributed uniform gravity grid database suitable for interpolation calculations. The indexing process includes the following steps:

[0071] Step 201: Select a corner point from the gravity database;

[0072] Assuming the selected gravity database has a latitude difference of 5′, the longitude difference is 10′ for areas above 60° north and south latitude, and 5′ for areas within 60° north and south latitude. On the three data blocks, select three points as corner points (e.g., ...). Figure 4 (Points shown within the dashed circle). The corresponding corner point information is: the number of points in this data block and the longitude interval.

[0073] Step 202: Calculate the required grid points using index interpolation;

[0074] If we interpolate a point on a plane that is not on a grid line (e.g.) Figure 4 The gravity value (shown by the gray dots) requires data from 4 grid points.

[0075] To find the data for these four grid points, we first need to determine the latitude block where the carrier is located. Assume... Figure 4 The coordinates of the gray point are [λ φ h] (where λ is latitude in degrees; φ is longitude in degrees; and h is altitude in meters). The steps are as follows:

[0076] 1) Calculate λ minus the latitude of corner point 2, that is: calculate λ-60 degrees;

[0077] 2) If λ-60 degrees > 0, then calculate λ minus the latitude of corner point 1, that is: calculate λ-90 degrees;

[0078] 3) If λ-90 degrees < 0, then the point is located on a latitude block within the range of 60 to 90 degrees;

[0079] Secondly, it is necessary to find the four grid points required for interpolation within the latitude block ranging from 60 to 90 degrees. The steps are as follows:

[0080] 1) Calculate 90 degrees - λ;

[0081] 2) Calculate the latitude interval Δλ, Δλ = floor((90 degrees - λ) × 60 / 5′), where floor represents rounding down;

[0082] 3) Calculate the longitude interval Δφ, Δφ=floor((φ-longitude of corner point 1)×60×cos60° / 5′;

[0083] 4) In a computer, if the data block corresponding to corner point 1 is stored as a one-dimensional vector (starting from the corner point 0, and stored sequentially from top to bottom and then from left to right), then the index point position of the interpolation point is calculated as follows:

[0084] point1 = 0 + Δφ × N + Δλ

[0085] point2 = 0 + Δφ × N + Δλ + 1

[0086] point3 = 0 + (Δφ + 1) × N + Δλ

[0087] point4 = 0 + (Δφ + 1) × N + Δλ + 1

[0088] Where N is (90 degrees - 60 degrees) × 60 / 5′.

[0089] Example 1

[0090] The existing gravity database interpolation scheme requires 106MB of storage space, which is relatively large and can be further optimized. The existing gravity grid stores data points in a uniform latitude and longitude format, with the density increasing closer to the poles. This approach has two drawbacks: 1. The data point density is high in high-latitude regions (especially near the poles), consuming a significant amount of storage space; 2. This storage method does not satisfy the fundamental assumptions of the interpolation algorithm. The interpolation algorithm assumes that gravity anomalies and vertical deviations occur within... Figure 5 The grid lines shown are distributed in a regular pattern. However, when a uniform latitude and longitude grid is selected for interpolation calculations, the longitude grid lines will be distorted in high-latitude regions (e.g., ...). Figure 6 (As shown).

[0091] The spherical uniform grid gravity compensation scheme divides the Earth's surface into 11 latitudinal zones. Within each latitudinal zone, the latitudinal intervals are the same, but the longitude intervals are different (e.g., ...). Figure 7 (As shown). The input to this scheme is the longitude, latitude, and elevation information calculated by the inertial navigation system. The algorithm logic is as follows:

[0092] 2.1 Obtain the solution output of the inertial reference system.

[0093] The required outputs of the inertial reference system are: the longitude (Lon, in degrees), latitude (Lat, in degrees), and elevation information (h, in meters) of the carrier.

[0094] 2.2 Pre-store database corner points.

[0095] If δ is set to 5 arcminutes and N_1 is set to 5, then 11 data blocks will be formed. The latitude intervals are randomly selected as follows: 90.0, 86.0, 81.0, 76.0, 71.0, 60.0, -60.0, -71.0, -76.0, -81.0, -86.0, -90.0. The latitude information of the data blocks is: double CORNER_BLOCK

[12] ={90.0, 86.0, 81.0, 76.0, 71.0, 60.0, -60.0, -71.0, -76.0, -81.0, -86.0, -90.0}; the corner information is: the longitude difference of each data block (in arcminutes, see Data_Block). double Data_Block

[11] = {60.0, 30.0, 20.0, 15.0, 10.0, 5.0, 10.0, 15.0, 20.0, 30.0, 60.0}; data points int INI_FLAG_BLOCK

[11] = {0, 17689, 61670, 127611, 215512, 502925, 6729486, 7016899, 7104800, 7170741, 7214722};

[0096] 2.3 Search for the latitude zone where the carrier is located.

[0097] (1) Pre-store latitude angles for different latitude zones ( Figure 7 The database contains the top left corner of each latitude point and the memory storage location of each corner point, see CORNER_BLOCK

[12] and INI_FLAG_BLOCK

[11] respectively.

[0098] double CORNER_BLOCK

[12] ={90.0,86.0,81.0,76.0,71.0,60.0,-60.0,-71.0,-76.0,-81.0,-86.0,-90.0};

[0099] int INI_FLAG_BLOCK

[11] ={0,17689,61670,127611,215512,502925,6729486,7016899,7104800,7170741,7214722};

[0100] (2) Set the latitude range of the initial search region to CORNER_BLOCK[init_hih_lat_flag] to CORNER_BLOCK[init_low_lat_flag] (where init_hih_lat_flag=5, init_low_lat_flag=6).

[0101] (3) Calculate the values ​​of res1 = lat - CORNER_BLOCK[init_hih_lat_flag] and res2 = lat - CORNER_BLOCK[init_low_lat_flag].

[0102] (4) When res1*res2>0 and res1<0, adjust the search region and let:

[0103] init_hih_lat_flag=init_hih_lat_flag+1

[0104] init_low_lat_flag=init_low_lat_flag+1

[0105] When res1*res2>0 and res1>0, adjust the search region as follows:

[0106] init_hih_lat_flag=init_hih_lat_flag-1

[0107] init_low_lat_flag=init_low_lat_flag-1

[0108] (5) If res1*res2<0, exit the loop, which means that the corresponding latitude zone has been found.

[0109] 2.4 Finding the interpolation point.

[0110] like Figure 8 , Figure 9 As shown, the interpolation point positions are:

[0111] point_LU=cor_pos+(lon_rd)*((int)(Delta_lat_block)+1)+lat_ld;

[0112] point_LD=cor_pos+(lon_rd)*((int)(Delta_lat_block)+1)+(lat_ld+1);

[0113] point_RU=cor_pos+(lon_rh)*((int)(Delta_lat_block)+1)+lat_ld;

[0114] point_RD=cor_pos+(lon_rh)*((int)(Delta_lat_block)+1)+(lat_ld+1);

[0115] Calculate the latitude and longitude distances a and b:

[0116] b = res2 / (5.0 / 60.0);

[0117] b = b - (int)b;

[0118] a=abs(lon+180.0) / (Data_Block[init_hih_lat_flag] / 60.0);

[0119] a = a - (int)a;

[0120] in:

[0121] cor_pos=INI_FLAG_BLOCK[init_hih_lat_flag]

[0122] lon_rd = (int)(abs(lon+180.0) / (Data_Block[init_hih_lat_flag] / 60.0))(lon is not equal to 180.0)

[0123] lon_rd = (int)(abs(lon+180.0) / (Data_Block[init_hih_lat_flag] / 60.0))-1 (lon equals 180.0)

[0124] lon_rh = lon_rd + 1;

[0125] Delta_lat_block=CORNER_BLOCK[init_hih_lat_flag]-CORNER_BLOCK[init_low_lat_flag] / (5 / 60.0)

[0126] lat_ld = (int)(abs(res1) / (5 / 60.0))(lat ≠ -90.0)

[0127] lat_ld = (int)(abs(res1) / (5 / 60.0)) - 1 (lat equals -90.0)

[0128] 2.5 Calculate gravity interpolation.

[0129] According to the method described in Section 2.4, the corresponding interpolation points are found in the storage space and named as G_0[4][3] and G_10000[4][3], respectively. The gravity value at the location of the carrier can be obtained by using fifth-order linear interpolation.

[0130] (1) Overview of Gravity Database

[0131] The number of database points in the 11 global latitude zones are as follows:

[0132] Num_1 = 361 × 49 = 17689;

[0133] Num_2 = 721 × 61 = 43981;

[0134] Num_3 = 1081 × 61 = 65941;

[0135] Num_4 = 1441 × 61 = 87901;

[0136] Num_5 = 2161 × 133 = 287413;

[0137] Num_6=4321×1441=6226561;

[0138] Num_7 = 2161 × 133 = 287413;

[0139] Num_8 = 1441 × 61 = 87901;

[0140] Num_9 = 1081 × 61 = 65941;

[0141] Num_10 = 721 × 61 = 43981;

[0142] Num_11 = 361 × 49 = 17689;

[0143] Total points num_total = 7232411.

[0144] Each grid point stores six short integer data points, including the gravity value in the E direction at 0 meters, the gravity value in the N direction at 0 meters, the gravity difference in the U direction at 0 meters (the difference between EGM2008 and the normal gravity model), the gravity value in the E direction at 10000 meters, the gravity value in the N direction at 10000 meters, and the gravity difference in the U direction at 10000 meters (the difference between EGM2008 and the normal gravity model).

[0145] The database is then stored in the following format:

[0146]

[0147] Short int type database size = num_total × 3 × 2 × 2 = 86788932 bytes ≈ 83MB.

[0148] Database file name: database_5_normal_short_int_0_10000_ion_type.bin.

[0149] (2) Horizontal interpolation implementation method

[0150] LSB = 1.0E-6

[0151]

[0152]

[0153]

[0154] Similarly, the following results can be obtained:

[0155]

[0156]

[0157]

[0158]

[0159]

[0160]

[0161]

[0162]

[0163]

[0164]

[0165]

[0166]

[0167]

[0168]

[0169]

[0170]

[0171]

[0172]

[0173]

[0174]

[0175]

[0176] The horizontal interpolation result of the random interpolation point is:

[0177]

[0178]

[0179]

[0180]

[0181]

[0182]

[0183] λ is the latitude value of the random interpolation point, in degrees, ranging from -90 to 90; λ is the longitude value of the random interpolation point, in degrees, ranging from -180 to 180.

[0184] (3) Method for interpolation in elevation direction

[0185]

[0186]

[0187]

[0188] Where h is the height of the random interpolation point, in meters (m).

Claims

1. A method for generating and indexing a grid database for gravity compensation interpolation calculation, characterized in that, The process includes the following steps: S1: Generate global gravity grid data using a high-precision global gravity model; S2: Based on the Earth's spherical variations, trim the latitudinal directions of the high-latitude grid in the generated global gravity grid data to generate a uniformly distributed gravity grid; S3: Determine the least significant bit of the data according to the requirements of the inertial reference system, and store the uniformly distributed gravity grid data in binary form to form a uniform gravity grid database; S4: Divide the data into latitudinal data blocks and store the information of each data block; S5: Index the data block where the carrier is located based on the information of each data block and the carrier's position coordinates; S6: Calculate the interpolation grid point data required to determine the gravity value at the carrier's location on the data block where the carrier is located. In step S1, the process of generating global gravity grid data specifically includes: selecting longitude and latitude at fixed intervals. The latitude and longitude grid of the angle minute is used to calculate and record the gravity field data at the intersection of latitude and longitude lines using the global ultra-high-order Earth gravity field model EGM-2008. In step S2, latitude In that region, the global gravity grid was cropped, and the longitude intervals of the gravity database were adjusted. Adjusted to: In step 3, the least significant bit of data is determined according to the requirements of the inertial reference system. The requirements of the inertial reference system refer to the zero-bias repeatability of the accelerometer, and the least significant bit refers to the number of significant digits that need to be retained. In step S4, for latitude The latitude zone is symmetrically divided into 2*N_1 data blocks along the equatorial plane; and The latitude data blocks together form a total of (2*N_1)+1 data blocks; the data block information includes the upper latitude boundary, the lower latitude boundary, the number of grid points of the data block, and the longitude interval; In step S5, the coordinates of the carrier position are: ,in Latitude, unit is degree; Longitude, in degrees; Height, unit: m ; The upper and lower latitude boundaries of each data block information obtained in step S4 are denoted as the dataset CORNER_BLOCK[i]. k上 i k下 ], where i k上 Indicates the upper dimensional boundary of the k-th data block; i k下 The upper dimension represents the lower boundary of the k-th data block; k = 1, ..., 2*N_1+1; i k上 = i （k-1）下 ; i k下 = i （k+1）上 i 1上 =90; i （2 N_1+1）下 =-90; The number of grid points in each data block is denoted as DATA_num = [NN] k ], k=1, ..., 2*N_1+1; where, NN k This represents the number of grid points in the k-th data block; make ; like Then continue the calculation. ,like Then the carrier is determined to be on the data block; if Then continue to order until found ; like Then let Continue calculating ,like Then calculate ,like If so, then the carrier is determined to be on that database block; if Then continue to order until found Determine the data block number k where the carrier is located, denoted as ; In step S6, let the larger absolute value of the upper and lower latitude boundaries of the data block containing the carrier be CORNER_BLOCK_MAX, and the smaller absolute value be CORNER_BLOCK_MIN. Then calculate... , where | represents the absolute value of the data; Calculate latitudinal intervals , ,in Indicates rounding down; Calculate longitude interval , ; in, The starting longitude of the data block where the carrier is located; The interpolation grid points are calculated as follows: in for .

2. A computer-readable storage medium storing a computer program, characterized in that, When the computer program is executed by a processor, it implements the steps of the method as described in claim 1.

Images