A method for measuring vertical deviation based on path-measured gravity anomaly compensation

Abstract

The present invention discloses a vertical deviation measurement method based on path-measured gravity anomaly compensation. Based on the "removal recovery" theory, the gravity anomaly integral region is first gridded, and the gravity anomaly residual value is calculated by interpolation using the measured gravity anomaly value on the carrier path. Finally, the high-frequency component of the vertical deviation is solved by integrating the Vening formula, and the low-frequency component calculated by the model is restored. The method includes: constructing a spherical harmonic model of the gravity field, calculating the long-wave component of the regional gravity anomaly and the vertical deviation; measuring the gravity anomaly on the path, designing a path determination interpolation algorithm, and interpolating the gravity anomaly residual; discretizing the Vening‑Meinesz formula, calculating the high-frequency component of the vertical deviation residual, and restoring the three steps of the long-wave component of the vertical deviation. The present invention uses the measured gravity anomaly data compensation on the carrier trajectory to directly solve the problem of obtaining a large amount of high-precision gravity anomaly data, greatly improves the measurement efficiency, and ensures the measurement accuracy of the vertical deviation.

Description

Technical Field

[0001] The invention belongs to the fields of geophysics, inertial navigation, geodesy and precision instrument technology, and in particular relates to a vertical line deviation measurement method based on path measured gravity anomaly compensation. Background Art

[0002] Due to the uneven distribution of mass inside the earth and the undulating surface, there is a large difference between the real gravity field of the earth and the normal gravity field of the artificially established reference ellipsoid. This difference is called the gravity disturbance field. The vertical deviation refers to the angle at which the real gravity vector deviates from the normal gravity vector, that is, the angle deviation between the vertical line on the geoid and the normal line of the reference ellipsoid.

[0003] At present, the direct differential method using inertial sensors is one of the methods to solve the compensation of vertical deviation and gravity disturbance. There are two major problems in directly measuring the gravity vector using inertial sensors: one is to deduct the influence of interference acceleration; the other is to provide a high-precision horizontal reference. The inertial navigation system can directly sense the specific force information through the accelerometer. The specific force refers to the difference between the absolute acceleration of the carrier relative to the inertial space and the gravitational acceleration, and can autonomously realize attitude measurement. According to Newton's second law, the specific force measured by the inertial sensor accelerometer contains the carrier motion acceleration information and gravity acceleration information, and the two cannot rely on inertial sensors to realize the separation of gravity information and acceleration information. Therefore, the effective separation of the carrier motion acceleration from the specific force information is the key to gravity measurement. At present, the carrier phase differential GPS technology can be used to measure the carrier's motion speed and acceleration with high precision, eliminating the influence of interference acceleration on inertial measurement. On the other hand, since the attitude error of the inertial navigation system is affected by the inertial device error and gravity error, and accumulates over time, the use of GPS position and velocity information and inertial navigation system for combined navigation can realize the estimation and compensation of system errors and improve the accuracy of attitude measurement. The Chinese patents with application numbers CN201910258966.8 and CN201910000427.4 have designed a gravity gradient measurement device and a real-time calculation method for vertical deviation based on horizontal tensor gravity gradient. The Chinese patents with application numbers CN201310730904.5 and CN202111056554.X use the attitude differential information of inertial navigation and reference attitude to solve the vertical deviation based on the attitude differential method, and respectively design a dynamic measurement device based on astronomical attitude maintenance and a system and terminal for attitude differential compensation. The Chinese patent with application number CN201410305314.2 designs a dynamic measurement device combining INS and GPS. The papers “A Method of Deflection of the Vertical Measurementbased on Attitude Difference Compensation” and “An Improved Method for Dynamic Measurement of Deflections of the Vertical Based on the Maintenance of Attitude Reference” study the estimation and compensation problems of vertical deviation based on the attitude difference method.

[0004] In addition to the direct difference method and the attitude difference method, the astronomical geodetic method, the gravity measurement method and the optimal estimation method are also means to effectively estimate the vertical deviation. The astronomical geodetic method is based on the observation of astronomical information to solve the astronomical longitude and latitude, and the vertical deviation is solved geometrically with the geodetic longitude and latitude. It has high accuracy, but the measurement method is single, time-consuming, and the measurement efficiency is not high. The Chinese patents with application numbers CN201510259859.9 and CN201710917104.2 have studied the astronomical geodetic method. The optimal estimation method is to statistically model the vertical deviation. Usually, the vertical deviation on the path is established as a first-order or second-order Markov model, and the gravity disturbance or vertical deviation is used as the system state quantity for optimal estimation. The papers “Research ongravity vertical deflection on attitude of position and orientation systemand compensation method” and “Dynamic measurement of high-frequency deflectionsof the vertical based on the observation of INS / GNSS integration attitudeerror” respectively establish the gravity disturbance information as first-order and second-order Markov models for filtering estimation.

[0005] Gravity measurement mainly uses the Vening formula to integrate the gravity anomaly information in the region on the spherical surface to solve the vertical deviation of the calculation point. In 1928, Vening and Meinesz derived the spherical integral expression for calculating the vertical deviation using regional gravity anomaly data, which is called the Vening-Meinesz formula. It is the basic formula for calculating the vertical deviation. At present, many scholars at home and abroad have improved this formula. The papers "Spherical Convolution Form of Vening-Meinesz Formula", "Inversion of Global Ocean Gravity Anomaly Based on Laplace Equation Vertical Deviation Method" and "The modified integral method for the determination of gravity disturbance near the Earth's surface" studied the calculation method of vertical deviation and gravity anomaly inversion based on Vening-Meinesz formula and inverse Vening-Meinesz formula. A technical difficulty in solving the vertical deviation using gravity anomaly data is the need to obtain a large amount of high-quality gravity anomaly data, which puts high demands on measurement efficiency and accuracy.

[0006] In summary, the current measurement methods of vertical deviation mainly include the following problems and technical defects:

[0007] 1. The direct difference method relies heavily on GNSS satellite navigation information, which requires high-precision acceleration information. However, satellite navigation signals are very susceptible to interference, and GNSS signals may be lost in tunnels and mountainous areas.

[0008] 2. The attitude differential method requires high-precision attitude reference information, which places high demands on the accuracy and stability of the gyroscope, and the accuracy of the vertical deviation measurement is easily disturbed.

[0009] 3. Gravity measurement requires a large amount of high-precision gravity anomaly data. The current single-point gravity anomaly mapping has a large workload and low efficiency, and it is difficult to map high-precision data over a large area. Summary of the invention

[0010] In view of the above-mentioned technical status, the purpose of the present invention is to provide a vertical line deviation measurement method based on path measured gravity anomaly compensation, so as to solve the problem of purely autonomous solution of vertical line deviation without relying on GNSS satellite guidance information and the difficulty of obtaining a large amount of high-precision gravity anomaly data by gravity measurement method.

[0011] The present invention and technical solutions are described as follows:

[0012] The basic concept of the present invention is to first grid the gravity anomaly integration area, and use the high-precision, high-resolution gravity field spherical harmonic model EIGEN-6C4 to solve the low-frequency information of gravity anomaly at the center point of each grid and the low-frequency component of the vertical deviation at the calculation point. Then, the measured gravity anomaly on the carrier path is used to compensate the gravity anomaly database value in the area, and finally, the high-frequency component of the vertical deviation is solved by the Vening formula integration, and the low-frequency component solved by the model is restored.

[0013] According to the above invention purpose and concept, the present invention provides a vertical deviation measurement method based on path measured gravity anomaly compensation, characterized in that: the gravity anomaly and the vertical deviation are respectively decomposed into the long-wave component of the earth's gravity field model and the medium- and short-wave components of the ground data for solution, the long-wave component of the model is removed by using the gravity anomaly measured data on the carrier trajectory, the remaining gravity anomaly residual, that is, the high-frequency component is interpolated, and integrated on the spherical surface to solve the residual amount of the vertical deviation, that is, the high-frequency component of the vertical deviation, and the long-wave component of the vertical deviation calculated by the gravity field model is restored at the calculation point, and finally the vertical deviation measurement method based on path measured gravity anomaly compensation is formed, comprising the following steps:

[0014] Step 1: Construct a spherical harmonic model of the gravity field and calculate the long-wave components of the regional gravity anomaly and vertical deviation;

[0015] Step 2: Measure gravity anomalies along the path, design a path determination interpolation algorithm, and interpolate gravity anomaly residuals;

[0016] Step 3: Discretize the Vening-Meinesz formula, calculate the high-frequency component of the vertical deviation residual, and restore the long-wave component of the vertical deviation;

[0017] The present invention further provides a method for measuring vertical deviation based on path measured gravity anomaly compensation, characterized in that: the specific steps of "constructing a spherical harmonic model of the gravity field, calculating the regional gravity anomaly and the long-wave component of the vertical deviation" in step 1 are as follows:

[0018] Step 1.1: Select the vertical deviation calculation point, determine the integration area, and perform grid processing;

[0019] Step 1.2: Construct the EIGEN-6C4 high-precision gravity field spherical harmonic model and calculate the gravity anomaly value Δg at the center point of each grid area EIG , and enter the latitude of the vertical deviation solution point Longitude λ c and height h c , calculate the long wave component of the vertical deviation of the solution point and the meridian component η EIG ,ξ EIG ;

[0020] Step 1.3: Perform initial alignment of the stationary base and input the coordinates of the carrier starting point, latitude Longitude λ, altitude h.

[0021] The present invention further provides a method for measuring vertical deviation based on path measured gravity anomaly compensation, characterized in that: the specific steps of "path measured gravity anomaly, designing path determination interpolation algorithm, and interpolating gravity anomaly residual" in step 2 are as follows:

[0022] Step 2.1: After the initial alignment is completed, the carrier enters the navigation mode, and the gravity measurement equipment performs dynamic gravity anomaly measurement calculations on the points on the carrier's trajectory and outputs the longitude and latitude of the path measurement points. and gravity anomaly Δg m ;

[0023] Step 2.2: Schematic diagram of residual interpolation algorithm is as follows Figure 2 As shown, it is determined whether the measured path passes through a single grid area. If it passes through a single grid area, steps 2-3 are executed; if it does not pass through the grid area, as shown in FIG. Figure 3 Grid area 3, then execute steps 2-5. The residual interpolation diagram is as follows Figure 3, as shown in Figure 4, the black center point in the figure is the center point of the integral grid area, the black L curve is the measured path, the black circle is the threshold area Ω, the point P is the vertical deviation calculation point, and the shaded grid area is the grid area passed by the measured path;

[0024] Step 2.3: Determine whether the current measured point crosses the threshold area Ω. If it crosses the threshold area, Figure 3 Grid area 1, the gravity anomaly measurements on the measured paths within the threshold area are averaged, and the average value on the path within the threshold is used as the measured gravity anomaly value of the center point in the grid area; if it does not pass through the threshold area, such as Figure 3 For grid area 2, Gaussian interpolation is performed on the gravity anomaly values ​​on the measured path within the grid area to obtain the measured gravity anomaly value at the center point of the grid area;

[0025] Step 2.4: Remove the long-wave component of the gravity anomaly calculated by the model from the measured value of the gravity anomaly at the center point of the grid area to solve the gravity anomaly residual;

[0026]

[0027] In the formula, represents the model-calculated gravity anomaly at the center point of the grid area traversed by the measured path; Represents the high-frequency component of the gravity anomaly residual at the center point.

[0028] Step 2.5: For the grid areas that are not passed by the entire measured path, such as Figure 4 In the white grid area, inverse distance weighted interpolation is performed to calculate the gravity anomaly residuals at the center points of the remaining grids.

[0029] The present invention further provides a method for measuring vertical deviation based on path measured gravity anomaly compensation, characterized in that: the specific steps of "discretizing the Vening-Meinesz formula, calculating the high-frequency component of the vertical deviation residual, and restoring the long-wave component of the vertical deviation" in step 3 are as follows:

[0030] Step 3.1: Calculate the Stokes function S(ψ), the spherical angular distance ψ, and the azimuth angle α respectively;

[0031]

[0032]

[0033]

[0034] In the formula, λ′ represents the latitude and longitude of the points in the integration area; λ represents the latitude and longitude of the vertical deviation calculation point.

[0035] Step 3.2: Schematic diagram of spherical integral is as follows Figure 5 As shown, the Vening-Meinesz formula (5) is discretized, that is, the double integral of the spherical surface is converted into the digital integral form (6), where the differential region dσ can be discretized into Δσ. Formula (6) can be further discretized into formulas (8) and (9). The integrals of each grid are summed to calculate the residual value ξ of the vertical deviation at the calculation point. res , η res ;

[0036]

[0037]

[0038]

[0039]

[0040]

[0041] In the formula, Represents the gravity anomaly residual value at the center point of the i-th grid area.

[0042] Step 3.3: Combine the vertical deviation residual value calculated in step 3.2 with the vertical deviation long-wave component calculated by the gravity field model to restore the vertical deviation long-wave component and finally solve the vertical deviation value at the calculation point.

[0043] ξ=ξ res +ξ EIG (10)

[0044] η=η res +η EIG (11)

[0045] The advantages of the present invention over the prior art are: without relying on GNSS satellite navigation information, the invention adopts an autonomous and efficient method for solving the vertical deviation, which solves the problem of directly obtaining a large amount of high-precision gravity anomaly data, greatly improves the measurement efficiency, and ensures the measurement accuracy of the vertical deviation. Compared with the traditional Vening-Meinesz formula based on the global area integration method, the accuracy of this method is greatly improved. BRIEF DESCRIPTION OF THE DRAWINGS

[0046] Figure 1 :Flowchart of a vertical line deviation measurement method based on path measured gravity anomaly compensation

[0047] Figure 2 :Schematic diagram of gravity anomaly residual interpolation algorithm

[0048] Figure 3 : Path interpolation diagram

[0049] Figure 4 : Path determination diagram

[0050] Figure 5 : Schematic diagram of spherical surface integral DETAILED DESCRIPTION

[0051] The specific implementation of the method of the present invention is now further described in detail with reference to the accompanying drawings, appendices and embodiments.

[0052] Example

[0053] Several simulation tests were conducted in Xi'an, Shaanxi Province (N 33°~N 36°, E 107.5°~110.5°), and several shipborne gravity anomaly actual tests were conducted in a certain area of ​​the South China Sea. The gravity measurement system was installed in the cabin, as shown in the attached picture. The vertical deviation measurement method based on path measured gravity anomaly compensation proposed by the present invention was experimentally verified. The technical solution of the present invention is described in detail by embodiments in conjunction with the accompanying drawings.

[0054] Step 1: Construct a spherical harmonic model of the gravity field and calculate the long-wave components of the regional gravity anomaly and vertical deviation;

[0055] Step 1.1: Determine a 3°×3° rectangular integration area (N 33°~N 36°, E 107.5°~110.5°) centered on the vertical deviation calculation point (N 34.5°, E 109°), and grid the area with a grid size of 2.5′×2.5′;

[0056] Step 1.2: Construct a high-precision, high-resolution EIGEN-6C4 gravity field spherical harmonic model. Use the model to calculate the gravity anomaly data Δg at the center point of each grid area. EIG , i.e. low frequency component, and input the latitude of the vertical deviation solution point (N34.5°), longitude λ c (E 109°) and elevation h c (430m), calculate the meridian component η of the long wave component of the vertical deviation at the solution point EIG ,ξ EIG ;

[0057] Step 1.3: Perform initial alignment of the stationary base and input the coordinates of the carrier starting point, latitude (N 33°25′47″), longitude λ(E 107°54′14″), altitude h(410m).

[0058] Step 2: Measure gravity anomaly along the path and interpolate gravity anomaly residuals;

[0059] Step 2.1: After the initial alignment is completed, the carrier enters the navigation mode, and the gravity measurement equipment performs dynamic gravity anomaly measurement calculations on the points on the carrier's trajectory and outputs the longitude and latitude of the path measurement points. and gravity anomaly Δg m ;

[0060] Step 2.2: Calculate the gravity anomaly residual value according to steps 2.2 to 2.5 in the invention content

[0061] Step 3: Calculate the vertical line deviation residual value (high-frequency component) based on the Vening-Meinesz formula, and restore the vertical line deviation long-wave component;

[0062] Step 3.1: Calculate the Stokes function S(ψ), the spherical angular distance ψ, and the azimuth angle α according to equations (2) to (4);

[0063] Step 3.2: Discretize the Vening-Meinesz formula according to equations (5) to (9), and sum the integrals of each grid to calculate the vertical line deviation residual value ξ at the calculation point res , η res ;

[0064] Step 3.3: Substitute the vertical line deviation residual value η calculated in step 3.2 res ,ξ res The long-wave component η of the vertical deviation calculated from the gravity field model EIG ,ξ EIG Combined, as in equations (10) and (11), the long-wave component of the vertical deviation is restored, and finally the vertical deviation values ​​η and ξ at the calculation point are solved;

[0065] Table 1 Calculated values ​​of vertical deviation at calculation points in simulation area

[0066]

[0067]

[0068] Table 2 Calculated values ​​of vertical deviation at calculation points in the test area

[0069]

[0070] The method for solving the vertical deviation proposed in the present invention is applied to simulation comparative analysis and actual experiments. The results of multiple simulations show that, in the selected area, the standard deviations of the meridian component and the vertical deviation are 1.042″ and 0.655″, respectively, and the root mean square errors are 1.679″ and 0.897″, as shown in Table 1. Compared with the traditional Vening-Meinesz formula based on integration in a global area, the accuracy of this method is greatly improved. The results of multiple tests in a certain sea area show that the standard deviations of the vertical deviation components solved by this method are 0.379″ and 1.437″, respectively, as shown in Table 2, which verifies the feasibility and accuracy of the method.

Claims

1. A method for measuring vertical deviation based on path measured gravity anomaly compensation, characterized in that: The gravity anomaly and the vertical deviation are respectively decomposed into the long-wave component of the earth's gravity field model and the medium- and short-wave components of the ground data for solution. The long-wave component of the model is removed by using the measured gravity anomaly data on the carrier trajectory. The remaining gravity anomaly residual, i.e., the high-frequency component, is interpolated and integrated on the spherical surface to solve the residual of the vertical deviation, i.e., the high-frequency component of the vertical deviation. The long-wave component of the vertical deviation calculated by the gravity field model is restored at the calculation point. Finally, a vertical deviation measurement method based on path measured gravity anomaly compensation is formed, which includes the following steps: Step 1: Construct a spherical harmonic model of the gravity field and calculate the long-wave components of the regional gravity anomaly and vertical deviation. The specific steps are as follows: Step 1.1: Select the vertical deviation calculation point, determine the integration area, and perform grid processing; Step 1.2: Construct the EIGEN-6C4 high-precision gravity field spherical harmonic model and calculate the gravity anomaly value Δg at the center point of each grid area EIG , and enter the latitude of the vertical deviation solution point Longitude λ c and height h c , calculate the long wave component of the vertical deviation of the solution point and the meridian component η EIG ,ξ EIG ; Step 1.3: Perform initial alignment of the stationary base and input the coordinates of the carrier starting point, latitude Longitude λ, altitude h; Step 2: Measure gravity anomalies along the path, design a path determination interpolation algorithm, and interpolate gravity anomaly residuals; Step 3: Discretize the Vening-Meinesz formula, calculate the high-frequency component of the vertical line deviation residual, and restore the long-wave component of the vertical line deviation.

2. A method for measuring vertical deviation based on path measured gravity anomaly compensation according to claim 1, characterized in that: The specific steps of "path-measured gravity anomaly, designing path determination interpolation algorithm, and interpolating gravity anomaly residual" described in step 2 are as follows: Step 2.1: After the initial alignment is completed, the carrier enters the navigation mode, and the gravity measurement equipment performs dynamic gravity anomaly measurement calculations on the points on the carrier's trajectory and outputs the longitude and latitude of the path measurement points. and gravity anomaly Δg m ; Step 2.2: According to the residual interpolation algorithm, determine whether the measured path passes through a single grid area. If it passes through a single grid area, execute step 2.3; if it does not pass through the grid area, execute step 2.5; Step 2.3: Determine whether the current measured point passes through the threshold area Ω. If it passes through the threshold area, then average the gravity anomaly measurement values ​​on the measured path within the threshold area, and use the average value on the path within the threshold as the measured gravity anomaly value of the center point of the grid area; if it does not pass through the threshold area, then Gaussian interpolation is performed on the gravity anomaly values ​​on the measured path within the grid area, so as to solve the gravity anomaly measured value at the center point of the grid area; Step 2.4: Remove the long-wave component of the gravity anomaly calculated by the model from the measured value of the gravity anomaly at the center point of the grid area to solve the gravity anomaly residual; In the formula, represents the model-calculated gravity anomaly at the center point of the grid area traversed by the measured path; Represents the high-frequency component of the gravity anomaly residual at the central point; Step 2.5: Perform inverse distance weighted interpolation on the grid areas that are not passed by the entire measured path to calculate the gravity anomaly residuals at the remaining grid center points.

3. The method for measuring vertical deviation based on path measured gravity anomaly compensation according to claim 1, characterized in that: The specific steps of "discretizing the Vening-Meinesz formula, calculating the high-frequency component of the vertical deviation residual, and restoring the long-wave component of the vertical deviation" described in step 3 are as follows: Step 3.1: Calculate the Stokes function S(ψ), the spherical angular distance ψ, and the azimuth angle α respectively; In the formula, λ' represents the latitude and longitude of the points in the integration area; λ represents the latitude and longitude of the vertical deviation calculation point; Step 3.2: According to the spherical integral formula, the Vening-Meinesz formula (5) is discretized, that is, the double integral of the spherical surface is converted into the digital integral form (6), where the differential region dσ can be discretized into Δσ. Formula (6) can be further discretized into formulas (8) and (9). The integrals of each grid are summed to calculate the residual value ξ of the vertical deviation at the calculation point. res , η res ; In the formula, represents the gravity anomaly residual value at the center point of the i-th grid area; Step 3.3: Combine the vertical deviation residual value calculated in step 3.2 with the vertical deviation long-wave component calculated by the gravity field model to restore the vertical deviation long-wave component, and finally solve the vertical deviation value at the calculation point. ξ=ξ res +ξ EIG (10) the=the res +n EIG (11).

Images