CN115165019A - Flight navigation satellite system (GNSS) tide level measuring method - Google Patents

Abstract

The invention discloses a method for measuring the sea-going GNSS sea level, which belongs to the technical field of navigation and is used for acquiring a sea level sequence and comprises the following steps: installing a measuring device, calculating a sea surface geodetic elevation sequence, extracting a tide level sequence based on geodetic elevation by adopting a variational modal decomposition and wavelet transformation combined denoising method, calculating the geodetic elevation of a seabed grid point, acquiring elevation information of a GPS leveling joint point in a measuring area, solving the approximate elevation abnormality of the seabed grid point by adopting a quadratic polynomial surface fitting method, summing the geodetic elevation and the approximate elevation abnormality to obtain the approximate normal height of the seabed grid point, solving the accurate elevation abnormality of the seabed grid point, solving the elevation abnormality of a measuring ship by adopting a bilinear interpolation method, and finally obtaining the tide level sequence based on the 1985 national elevation standard. Compared with the prior art, the tidal level root mean square error extracted by the noise elimination method provided by the invention is smaller and is closer to a real tidal level value. The coordinates and the ground height of the sea bottom point are accurately measured by using the sonar, and the accuracy of vertical reference conversion is improved.

Description

Flight navigation satellite system (GNSS) tide level measuring method

Technical Field

The invention discloses a method for measuring the sea level of an underway GNSS, and belongs to the technical field of navigation.

Background

In the process of measuring the submarine topography by using the multi-beam sonar, tide level data of an investigation area are synchronously obtained, and the step of completing tide level correction is an indispensable step of data processing. The traditional tide level measuring mode comprises a water gauge tide test, a digital pressure tide test instrument tide test and a long-term tide test station tide test. The water gauge tide checking method comprises the steps of vertically fixing the water gauge on a wharf wall and a beach, manually reading the position of a water surface on the water gauge, recording the height of a sea surface on the zero point of the water gauge, and timing by combining a clock; the digital pressure tide gauge tests tide by measuring the pressure change of seawater to indirectly calculate the lifting change of the sea surface; the long-term tide gauge station is used for measuring the tide level at a long-term fixed tide gauge point set in the country by using a water gauge or a tide gauge.

The traditional tide gauge mode is affected by many adverse factors such as tide level model error, no tide gauge station, loss of a tide gauge, consumption of a large amount of human resources, time consumption and the like. With the continuous improvement of the GNSS carrier phase differential measurement technology, the GNSS-RTK, GNSS-PPK and network RTK tide level measurement methods are gradually popularized and applied in China. The GNSS tide gauge is used for obtaining high-precision elevation information of a shipborne antenna by utilizing a GPS carrier phase technology, obtaining an elevation sequence of a sea surface through the height difference between the antenna and the sea surface, extracting low-frequency information in the elevation sequence of the sea surface by utilizing a filtering method, and obtaining the GNSS flight altitude through vertical reference conversion.

The existing GNSS tide gauge technical scheme comprises the steps of coordinate system establishment, GNSS receiver geodetic height data acquisition, sea surface geodetic height sequence filtering and vertical reference conversion, and the specific steps are as follows:

(1) Before measurement begins, a ship body coordinate system is established by taking the IMU as a reference origin, the ship bow direction is an X axis, the vertical X axis points to a starboard and is a Y axis, and the vertical XOY plane points downwards and is a Z axis. And expressing the coordinates of the two GNSS receivers in a ship body coordinate system.

(2) The longitude, the latitude and the geodetic elevation of a GNSS receiver are obtained by utilizing a GNSS-RTK or GNSS-PPK or network RTK technology, a geodetic altitude sequence of the sea surface is obtained by combining the coordinate of the GNSS receiver under the established coordinate system, the attitude of a ship body measured by an IMU and the dynamic draft measured by a pressure depth measuring sensor in the step (1), high-frequency noise in the geodetic altitude sequence of the sea surface is removed by a fast Fourier transform (a flow chart of a Fourier transform method is shown in figure 1) or a wavelet forced noise elimination method, and a tide level value under 1985 national elevation reference is obtained by means of elevation reference conversion, wherein the traditional mode of the elevation reference conversion is a geometric method or a parameter conversion method.

At present, the coordinate, the geodetic height and the average geodetic height of a seabed grid point measured by multi-beam precision are not considered in the conventional elevation reference conversion mode to further improve the precision of vertical reference conversion, a fast Fourier transform filtering method adopted by GNSS tide observation cannot process nonstationary signals and abrupt signal wavelet transform to force selection of noise elimination wavelet bases to be difficult and lack of adaptivity, and the filtering effects of the two filtering methods are poor. When the vertical reference conversion of the sea level is carried out in a complicated submarine topography area, the conversion precision of a single geometric method or a conversion parameter method is low.

Disclosure of Invention

The invention discloses a method for measuring the sea-going GNSS tide level, which aims to solve the problems of poor filtering effect and low conversion precision of a GNSS tide gauge method in the prior art.

A method for measuring an aviation GNSS sea level comprises the following steps:

s1, mounting a measuring device on a ship board, wherein the measuring device comprises a multi-beam sonar, an IMU (inertial measurement unit), a pressure depth sounding sensor and two GNSS (global navigation satellite system) receivers, and the IMU is an inertial measurement unit;

the method comprises the steps that a GNSS receiver independently receives satellite signals, the longitude and latitude and the geodetic elevation of a phase center of the GNSS receiver are obtained, sonar measures the underwater topography of a survey area, the longitude and latitude, geodetic height data and depth of an underwater survey point are obtained, an IMU obtains a heave value at the IMU and an attitude value of a survey ship, the attitude value comprises a roll value, a pitch value and a yaw value, and a pressure depth measurement sensor measures the distance from the IMU to the water surface in a dynamic draft state when the survey ship moves and the distance from the IMU to the water surface when the ship is anchored;

s2, calculating a sea surface geodetic elevation sequence under the WGS84 elevation reference by integrating the data obtained in the S1;

s3, extracting a sea surface geodetic elevation sequence based on a geodetic elevation by adopting a variational modal decomposition and wavelet transform combined denoising method;

s4, solving the geodetic elevation of the seabed grid points of the measuring area based on the WGS84 elevation standard;

s5, obtaining elevation information of 6 or more than 6 GPS leveling combined measuring points near a measuring area, wherein the elevation information simultaneously comprises geodetic elevation of the GPS leveling combined measuring points and elevation under 1985 national elevation standard;

s6, solving the approximate normal height of the seabed grid points in the measuring area by adopting a quadratic polynomial surface fitting method;

s7, solving the precise elevation abnormity of the seabed grid points according to the MoroKingki principle;

and S8, solving out the precise elevation abnormity at the position of the ship to be tested by adopting a bilinear interpolation method according to the longitude and latitude of the phase center of the GNSS receiver in the S1, and reducing the tide level sequence in the S3 to the 1985 national elevation standard to obtain the tide level sequence based on the 1985 national elevation standard.

Preferably, S1 comprises the following sub-steps:

s1.1, establishing an ideal hull coordinate system by taking an IMU as a reference origin, establishing a three-dimensional model of a measuring device in mapping software, measuring the offset relation between each component on the measuring device through the three-dimensional model, expressing the coordinates of two GNSS receivers under the ideal hull coordinate system, measuring the vertical distance from the reference origin to the GNSS receivers through the three-dimensional model and recording the vertical distance as L1, and obtaining the distance L2 from the IMU to the water surface when a hull is anchored through a pressure depth measuring sensor, namely obtaining the distance L3 from the GNSS receivers to the water surface, wherein the calculation formula is as follows: l3= L1-L2;

s1.2, during sonar measurement, obtaining a roll mounting deviation angle and a pitch mounting deviation angle of the transducer by adopting a patch test in a typical area;

s1.3, correcting the roll value and the pitch value measured by the IMU according to the roll installation deviation angle and the pitch installation deviation angle obtained in the S1.2, and correcting the pitch value measured by the IMU only if the sonar system has a roll compensation function;

and S1.4, acquiring accurate geodetic height and longitude and latitude data of two GNSS receivers in phase centers in a GNSS-PPK, GNSS-RTK, network RTK or PPP mode.

Preferably, S2 comprises the following sub-steps:

s2.1, respectively calculating induced heave values H of the two GNSS receivers according to the coordinates of the two GNSS receivers in the S1.1 in the ideal hull coordinate system and the pitch values and the roll values corrected in the S1.3 _i ：H _i ＝xsin P-ysin Rcos P-z(cos Rcos P-1)，H _i Representing induced heave values, x, y and z representing coordinates of the GNSS receiver in an ideal hull coordinate system, and P and R representing pitch values and roll values recorded by the IMU in S1;

s2.2, according to the induced heave value H in S2.1 _i And (5) calculating the heave value h of the two GNSS receivers respectively according to the heave value of the IMU in the S1 _G ：h _G ＝H _i +h ₀ ，h ₀ Represents the heave value at the IMU in S1;

s2.3, the calculation process of the sea surface geodetic elevation sequence is as follows: h = h _GPS -h _G -H _m + delta H, H is the instantaneous ground height of the sea surface, namely the sea surface ground height sequence; h is _GPS Geodetic elevation of a phase center obtained by the GNSS receiver in S1; h _m Is L3; Δ H is the difference between the dynamic draft of the pressure depth sensor measuring the vessel motion in S1 and L3 in S1.1.

Preferably, S3 comprises the following sub-steps:

s3.1, successively carrying out K-layer modal decomposition on the sea surface geodetic sequence obtained in the S2 by adopting a variational modal decomposition method, wherein K =10n, n =1,2,3 \ 8230, n is a positive integer, a penalty coefficient alpha =80000, and a convergence tolerance tol =1e-7;

s3.2, calculating signal energy E obtained by performing K-layer modal decomposition on the sea surface geodetic high sequence each time:

n is the number of sampling points, IMF _K (t) a K-th layer signal sequence for each modal decomposition;

s3.3, solving an energy ratio:

when the energy ratio gamma tends to be stable, the gamma fluctuates at the value of 0, and K is an appropriate modal decomposition number;

s3.4, decomposing the mode with the number of K to obtain a first layer mode IMF ₁ And performing variation modal decomposition again, wherein the penalty coefficient alpha =80000, the convergence tolerance tol =1e-7 and is called iterative decomposition for 1 time, the judgment criterion of the K value of the iterative decomposition for 1 time refers to S3.2 and S3.3, and the first layer of mode obtained by the iterative decomposition for 1 time is marked as IMF ₁₁ The second layer mode is denoted as IMF ₂₂ The IMF ₁₁ And IMF ₂₂ Reconstructing to obtain a new signal;

s3.5. For pair IMF in S3.4 ₁₁ And IMF ₂₂ Performing spectrum analysis to determine IMF ₁₁ And IMF ₂₂ The signal spectrum composition of (a);

s3.6. The maximum frequency omega corresponding to the sea surface earth elevation sequence in S2 is as follows:

f _s is the sampling frequency of the GNSS signal;

assuming that the frequency of other elements of the sea surface elevation sequence is f _w The frequency of the tide being f _t The number of layers z using the wavelet number of layers is determined by the following formula:

IMF of S3.4 ₁₁ And IMF ₂₂ Decomposing the z layer by the reconstructed new signal with sym4 wavelet to obtain the z-th layer low-frequency coefficient A _z And a high frequency coefficient D ₁ D ₂ ………D _z The frequency band occupied by the low frequency coefficient is the frequency band of the tideThe frequency band occupied by the high-frequency coefficient is the frequency band of other high-frequency elements, the high-frequency coefficient is set to zero and then is reconstructed with the low-frequency coefficient, so that the forced noise elimination of the sea surface elevation sequence can be realized, and the extracted low-frequency signal is a tide level sequence h based on WGS84 elevation reference _{tide_84} 。

Preferably, S4 comprises the following sub-steps:

s4.1, carrying out installation deviation angle correction, sound velocity correction and gross error elimination on the water depth data measured by sonar;

s4.2, adopting the h acquired in the S3 to the water depth data processed in the S4.1 _{tide_84} Carrying out tide level correction:

H _i ＝|h _i -h _{tide_84} |，H _i for the elevation of the subsea measurement point, h _i The water depth data after S4.1 processing is obtained;

s4.3, the sea surface geodetic elevation sequence in the S2 is subjected to grid formation according to 500m × 500m, and geodetic coordinates of grid points under a WGS84 geodetic coordinate system or plane coordinates and geodetic height h are output _Ω Weighting and averaging elevation values of grid points to obtain the average ground height under WGS84 elevation standard, and recording as h _γ 。

Preferably, S6 comprises the following sub-steps:

s6.1, assuming the existence of a quadratic polynomial surface fitting mathematical model between the elevation abnormity and the plane coordinates on the leveling joint measurement point:

a ₀ 、a ₁ 、a ₂ 、a ₃ 、a ₄ 、a ₅ for the model undetermined parameter, ζ _i When the number of the leveling joint measuring points is more than 6, the following error equation is formed for the elevation abnormity of the leveling joint measuring points: v = AX- ζ ⁰ The least square principle can be used to obtain A = (X) ^T X) ^-1 X ^T ζ ⁰

According to the minimumSolving the undetermined parameter a of the model by the solution of the two-times principle ₀ 、a ₁ 、a ₂ 、a ₃ 、a ₄ 、a ₅ The value of (d);

s6.2, substituting the coordinates of the seabed grid points in the S4 into the quadratic polynomial surface fitting mathematical model of the S6.1 to interpolate approximate elevation abnormity zeta of the seabed grid points ₀ ；

S6.3, measuring the approximate normal height h of the seabed grid points of the area _Γ Comprises the following steps: h is a total of _Γ ＝ζ ₀ +H _i ，ζ ₀ Indicating approximate elevation anomaly, H _i Is the geodetic elevation of the seabed grid point in S4.

Preferably, S7 comprises the following sub-steps:

s7.1 obtaining h obtained in S6 _Γ ；

S7.2, the elevation abnormity zeta of the seabed grid point is determined by a long wave part, namely the approximate elevation abnormity zeta ₀ And short wave part ζ _r Obtaining: ζ = ζ ₀ +ζ _r ；

S7.3, calculating short wave part zeta in elevation abnormity _r ：ζ _r ＝T _r /r，T _r Influence of topographic relief on disturbance positions of the submarine grid points; ζ represents a unit _r Correcting the terrain;

g is universal gravitation constant, G =6.673 × 10 ^-8 C ³ s ^-2 g ^-1 (ii) a ρ is the earth's mean mass density, ρ =2.67gC ^-3 ，γ ₀ ＝[(x-x _p ) ² +(y-y _p ) ² ] ^1/2 (x, y) are the coordinates of the sea-bottom grid point in S4, (x) _p ,y _p ) The coordinate of a certain grid point P in the seabed grid points is obtained;

s7.4. The normal gravity r similar to the ground level: r = r ₀ -0.3086h _Γ ，r ₀ For reference to normal gravity values on an ellipsoid, position: 10 ^-5 …m/s ² ，

Measuring point geographical latitude in a WGS84 geodetic coordinate system, unit: and (4) degree.

Preferably, S8 comprises the following sub-steps:

s8.1, according to the longitude and latitude of the phase center of the GNSS receiver in S1, the geodetic elevation of the seabed grid point in S4 and zeta in S7, adopting a bilinear interpolation method to solve zeta which is abnormal in elevation at the position of the ship to be measured ₁ ；

S8.2, calculating a finally interpolated elevation abnormal value f (x, y):

f(f,y ₁ )、f(x,y ₂ ) Showing the elevation anomaly, ζ, interpolated in the x-direction ₁₁ 、ζ ₂₁ 、ζ ₁₂ 、ζ ₂₂ Representing elevation outliers, x, at four corners of the grid ₁ 、x ₂ 、y ₁ 、y ₂ The coordinates of the four corner points in the X direction and the Y direction are represented, and the coordinates of the ship in the X direction and the Y direction are represented by X and Y;

s8.3, calculating a tide level sequence h based on 1985 national elevation standard _{tide_85} ：h _{tide_85} ＝h _{tide_84} +ζ ₁ 。

The invention has the main beneficial effects that: the signal-to-noise ratio of the sea level elevation sequence extracted by the new variational modal decomposition and wavelet transformation combined denoising method is high, namely the denoising method provided by the invention has good denoising effect, the root mean square error between the sea level and the real sea level extracted by the denoising method is small, namely the GNSS navigation tide level value extracted is closer to the real sea level value. The coordinate, the geodetic height and the average geodetic height of the sea bottom grid points accurately measured by the sonar can improve the precision of the vertical reference conversion.

Drawings

FIG. 1 is a flow chart of a Fourier transform method in the prior art;

FIG. 2 is a schematic structural diagram of a measuring apparatus;

FIG. 3 is a flow chart of extracting a sea level elevation sequence by using a variational modal decomposition and wavelet transform combined denoising method;

FIG. 4 is a flow chart of tide level reduction to 1985 national elevation reference based on WGS84 elevation reference;

FIG. 5 is a general flow chart of a GNSS-based method of performing an aviation sea level survey;

FIG. 6 is a schematic elevation view of a seafloor point based on WGS84 elevation reference;

FIG. 7 is a schematic diagram of calculating an elevation anomaly of a seafloor grid point;

FIG. 8 is a schematic illustration of calculating an elevation anomaly at a survey vessel;

FIG. 9 instantaneous ground height of the sea surface;

fig. 10 is a diagram illustrating extracting a tide bitmap in a sea surface elevation sequence based on a variational modal decomposition and wavelet transform joint denoising method in the embodiment of the present invention.

The reference numerals include: the device comprises a 1-GNSS receiver, a 2-receiver support rod, a 3-receiver support beam, a 4-receiver mounting plate, a 5-first upright post, a 6-upright post mounting plate, a 7-second upright post, a 8-third upright post, a 9-IMU, a 10-multi-beam transmitting transducer, a 11-multi-beam receiving transducer, a 12-multi-beam mounting plate and a 13-pressure sounding sensor.

Detailed Description

The following description will further illustrate embodiments of the present invention with reference to specific examples:

a method for measuring an aviation GNSS sea level comprises the following steps:

s1, mounting a measuring device on a ship board, wherein the measuring device comprises a multi-beam sonar, an IMU9, a pressure sounding sensor 13 and two GNSS receivers 1, and the IMU9 is an inertia measuring unit; IMU9 Inertial \8230, measurement \8230, unit, inertial Measurement Unit.

The method comprises the following steps that a GNSS receiver 1 independently receives satellite signals, the longitude and latitude and the geodetic elevation of a phase center of the GNSS receiver 1 are obtained, the sonar measures the underwater topography of a survey area, the longitude and latitude, geodetic height data and depth of an underwater survey point are obtained, an IMU9 obtains a heave value at the IMU9 and a ship survey attitude value, the attitude value comprises a roll value, a pitch value and a yaw value, and a pressure depth measurement sensor 13 measures the distance from the IMU9 to the water surface in a dynamic draft state when the ship surveys the motion and the distance from the IMU9 to the water surface when the ship is anchored;

s2, calculating a sea surface earth elevation sequence under the WGS84 elevation standard by integrating the data obtained in the S1;

s3, extracting a sea surface geodetic elevation sequence tidal level sequence based on geodetic elevation by adopting a variational modal decomposition and wavelet transform combined denoising method, as shown in FIG. 10;

s4, solving the geodetic elevation of the seabed grid points of the measuring area based on the WGS84 elevation standard;

s5, obtaining elevation information of 6 or more than 6 GPS leveling combined measuring points near a measuring area, wherein the elevation information simultaneously comprises geodetic elevation of the GPS leveling combined measuring points and elevation under 1985 national elevation standard;

s6, solving the approximate normal height of the seabed grid points in the measuring area by adopting a quadratic polynomial surface fitting method;

s7, solving the elevation abnormity of the seabed grid points according to the MoroKingki principle;

s8, solving out the accurate elevation abnormity at the ship measuring position by adopting a bilinear interpolation method according to the longitude and latitude of the phase center of the GNSS receiver 1 in the S1, and reducing the tide sequence in the S3 to the 1985 national elevation standard to obtain the tide sequence based on the 1985 national elevation standard.

S1 comprises the following substeps:

s1.1, establishing an ideal ship body coordinate system by taking the IMU9 as a reference origin, establishing a three-dimensional model of a measuring device in mapping software, measuring the offset relation among all parts on the measuring device through the three-dimensional model, expressing the coordinates of two GNSS receivers 1 under the ideal ship body coordinate system, measuring the vertical distance from the reference origin to the GNSS receivers 1 through the three-dimensional model and recording the vertical distance as L1, and obtaining the distance L2 from the IMU9 to the water surface when a ship body is anchored through a pressure depth measuring sensor 13, namely obtaining the distance L3 from the GNSS receivers 1 to the water surface, wherein the calculation formula is as follows: l3= L1-L2;

s1.2, during sonar measurement, obtaining a roll installation deviation angle and a pitch installation deviation angle of the transducer in a typical area by adopting a 'patch test';

s1.3, correcting the roll value and the pitch value measured by the IMU9 according to the roll installation deviation angle and the pitch installation deviation angle obtained in the S1.2, and correcting only the pitch value measured by the IMU9 if the sonar system has a roll compensation function;

and S1.4, acquiring accurate geodetic height and longitude and latitude data of the phase centers of the two GNSS receivers 1 in a GNSS-PPK, GNSS-RTK, network RTK or PPP mode.

S2 comprises the following substeps:

s2.1, respectively calculating induced heave values H of the two GNSS receivers 1 according to the coordinates of the two GNSS receivers 1 in the S1.1 in the ideal hull coordinate system and the pitch values and the roll values corrected in the S1.3 _i ：H _i ＝xsin P-ysin Rcos P-z(cos Rcos P-1)，H _i Representing induced heave values, x, y and z representing coordinates of the GNSS receiver 1 in an ideal hull coordinate system, and P and R representing pitch values and roll values recorded by the IMU9 in S1;

s2.2 according to the induced heave value H in S2.1 _i And the heave value at the IMU9 in S1, and respectively calculating the heave value h at the two GNSS receivers 1 _G ：h _G ＝H _i +h ₀ ，h ₀ Represents the heave value at IMU9 in S1;

s2.3, the calculation process of the sea surface earth elevation sequence is as follows: h = h _GPS -h _G -H _m + delta H, H is the instantaneous ground height of the sea surface, namely the sea surface ground height sequence; h is a total of _GPS Geodetic elevation of the phase center obtained by the GNSS receiver 1 in S1; h _m Is L3; Δ H is the difference between the dynamic draft measured by the tonometer sensor 13 in S1 during the survey of the ship and L2 of S1.1.

S3 comprises the following substeps:

s3.1, successively carrying out K-layer modal decomposition on the sea surface geodetic sequence obtained in the S2 by adopting a variational modal decomposition method, wherein K =10n, n =1,2,3 \ 8230, n and n are positive integers, the penalty coefficient alpha =5000, and the convergence tolerance tol =1e-7;

s3.2. Calculating K layers of sea surface geodetic height sequence every timeSignal energy E obtained by modal decomposition:

n is the number of sampling points, IMF _K (t) a K-th layer signal sequence for each modal decomposition;

s3.3, solving an energy ratio:

when the energy ratio gamma tends to be stable, the gamma fluctuates at the value of 0, and K is an appropriate modal decomposition number;

s3.4, decomposing when the mode decomposition number is K to obtain a first layer mode IMF ₁ And performing variation modal decomposition again, wherein the penalty coefficient alpha =80000, the convergence tolerance tol =1e-7 and is called iterative decomposition for 1 time, the K value judgment criterion of the iterative decomposition for 1 time refers to S3.2 and S3.3, and the first-layer mode obtained by the iterative decomposition for 1 time is recorded as IMF ₁₁ Second layer mode is denoted as IMF ₂₂ The IMF ₁₁ And IMF ₂₂ Reconstructing to obtain a new signal;

s3.5. For pairs IMF in S3.4 ₁₁ And IMF ₂₂ Performing spectrum analysis to determine IMF ₁₁ And IMF ₂₂ The signal spectrum composition of (a);

s3.6. The highest frequency omega corresponding to the sea surface earth elevation sequence in S2 is as follows:

f _s is the sampling frequency of the GNSS signal;

assuming that the frequency of other elements of the sea surface elevation sequence is f _w The frequency of the tide being f _t The number of layers z using the wavelet number of layers is determined by the following formula:

the IMF described in S3.4 ₁₁ And IMF ₂₂ Decomposing the z layer of the reconstructed signal by sym4 wavelet to obtain the low-frequency coefficient A of the z layer _z And a high frequency coefficient D ₁ D ₂ ………D _z The frequency band occupied by the low-frequency coefficient is the frequency band of tide, the frequency band occupied by the high-frequency coefficient is the frequency band of other high-frequency elements, the high-frequency coefficient is reset and then reconstructed with the low-frequency coefficient, the forced noise elimination of the sea surface elevation sequence can be realized, and the extracted low-frequency signal is the tide level sequence h based on the WGS84 elevation reference _{tide_84} 。

S4 comprises the following substeps:

s4.1, carrying out installation deviation angle correction, sound velocity correction and gross error elimination on the water depth data measured by sonar;

s4.2, adopting the h acquired in the S3 to the water depth data processed in the S4.1 _{tide_84} Carrying out tide level correction:

H _i ＝|h _i -h _{tide_84} |，H _i for the elevation of the subsea measurement point, h _i The water depth data after S4.1 processing is obtained;

s4.3, the sea surface geodetic elevation sequence in the S2 is subjected to grid formation according to 500m × 500m, and geodetic coordinates of grid points under a WGS84 geodetic coordinate system or plane coordinates and geodetic height h are output _Ω Weighting and averaging the elevation values of grid points to obtain the average ground height under the WGS84 elevation standard, and recording as h _γ 。

S6 comprises the following substeps:

s6.1, assuming the existence of a quadratic polynomial surface fitting mathematical model between the elevation abnormity and the plane coordinates on the leveling joint measurement point:

a ₀ 、a ₁ 、a ₂ 、a ₃ 、a ₄ 、a ₅ for the model to be determined parameter, ζ _i When the number of the leveling joint measuring points is more than 6, the following error equation is formed for the elevation abnormity of the leveling joint measuring points: v = AX- ζ ⁰ A = (X) is obtained according to the least square principle ^T X) ^-1 X ^T ζ ⁰

Solving the undetermined parameter a of the model according to least square primary understanding ₀ 、a ₁ 、a ₂ 、a ₃ 、a ₄ 、a ₅ The value of (d);

s6.2, substituting the coordinates of the seabed grid points in the S4 into the quadratic polynomial surface fitting mathematical model of the S6.1 to interpolate approximate elevation abnormity zeta of the seabed grid points ₀ ；

S6.3, measuring the approximate normal height h of the seabed grid points of the area _Γ Comprises the following steps: h is _Γ ＝ζ ₀ +H _i ，ζ ₀ Indicating an approximate elevation anomaly, H _i Is the geodetic elevation of the seabed grid points in S4.

S7 comprises the following substeps:

s7.1. Obtaining h obtained in S6 _Γ ；

S7.2, the elevation abnormity zeta of the seabed grid point is determined by a long wave part, namely the approximate elevation abnormity zeta ₀ And short wave part ζ _r Obtaining: ζ = ζ ₀ +ζ _r ；

S7.3, calculating short wave part zeta in elevation abnormity _r ：ζ _r ＝T _r /r，T _r Influence of topographic relief on disturbance positions of the submarine grid points is achieved; zeta _r Correcting the terrain;

g is universal gravitation constant, G =6.673 × 10 ^-8 C ³ s ^-2 g ^-1 (ii) a ρ is the earth's mean mass density, ρ =2.67gC ^-3 ，γ ₀ ＝[(x-x _p ) ² +(y-y _p ) ² ] ^1/2 (x, y) are the coordinates of the sea-bottom grid point in S4, (x) _p ,y _p ) The coordinate of a certain grid point P in the seabed grid points is obtained;

s7.4, normal gravity r similar to a ground level: r = r ₀ -0.3086h _Γ ，r ₀ For reference to normal gravity values on an ellipsoid, position: 10 ^-5 …m/s ² ，

Measuring point geographical latitude in a WGS84 geodetic coordinate system, unit: and (4) degree.

S8 comprises the following substeps:

s8.1, according to the longitude and latitude of the phase center of the GNSS receiver 1 in S1, the geodetic elevation of the seabed grid point in S4 and zeta in S7, adopting a bilinear interpolation method to solve zeta of elevation abnormity at the position of the ship to be measured ₁ ；

S8.2, calculating a finally interpolated elevation abnormal value f (x, y):

f(x,y ₁ )、f(x,y ₂ ) Showing the elevation anomaly, ζ, interpolated in the x-direction ₁₁ 、ζ ₂₁ 、ζ ₁₂ 、ζ ₂₂ Representing elevation outliers, x, at four corners of the grid ₁ 、x ₂ 、y ₁ 、y ₂ The coordinates of the four corner points in the X direction and the Y direction are represented, and the coordinates of the ship in the X direction and the Y direction are represented by X and Y;

s8.3, calculating a tide level sequence h based on 1985 national elevation standard _{tide_85} ：h _{tide_85} ＝h _{tide_84} +ζ ₁ 。

As shown in fig. 1, the measuring device is welded on a ship board of a ship survey through two mounting plates fixed on vertical columns, and the measuring device comprises two GNSS receivers 1, a receiver support rod 2, a receiver support beam 3, a receiver mounting plate 4, a first vertical column 5, a vertical column mounting plate 6, a second vertical column 7, a third vertical column 8, a multi-beam transmitting transducer 10, a multi-beam receiving transducer 11, an IMU9, a pressure depth sounding sensor 13 and a multi-beam mounting plate 12.

The GNSS receiver 1 receives signals transmitted by satellites using GNSS-PPK, GNSS-RTK, network RTK or

High-precision longitude and latitude and geodetic height data are obtained by a PPP mode, and the sampling frequency is 1 \ 8230HZ;

b. the multi-beam sounding sonar comprises a transmitting transducer and a receiving transducer, and obtains offshore high-precision submarine topography in a 'strip' covering measurement mode;

the IMU9 and the transducer are installed underwater, the IMU9 provides an attitude value and a heave value of the transducer, a ship body coordinate system is established by taking the IMU9 as a reference origin O, an X axis horizontally points forward to the bow direction (the measuring device P is parallel to the bow direction when being installed), a Y axis horizontally passes through the reference origin to the right, and a Z axis is vertical to the XOY plane and is vertically downward;

d. the sampling frequency of the pressure depth sounding sensor 13 is 5 \ 8230HZ, and the distance value from the IMU9 reference origin to the water surface is provided when the ship is anchored, so that the water depth measured by the transducer is corrected; providing dynamic draft of a survey ship when the ship moves so as to calculate the instantaneous ground height of the sea surface;

fig. 4 is a general flow chart of the method, fig. 2 is a flow chart of extracting a tide level (WGS 84 elevation reference) from the geodetic altitude of the sea surface by a variational modal decomposition and wavelet transform combined denoising method, and fig. 9 is a result chart thereof, and the number z of layers of wavelet decomposition is judged according to the tidal characteristics of a measurement area. Taking Qingdao semisolar tide sea area as an example, the tide cut-off frequency is about 1/(3600 x 3) ≈ 9.26 x 10 ^{- 5} Hz, the period of variation of high-frequency noise is a few seconds to a dozen minutes at most, and the lowest frequency of the noise is 1/1200 ≈ 8.33 × 10 ^-4 Hz, 0.5/2 when the number of layers z of the wavelet decomposition is taken to be 10 at minimum, in order to separate the tidal signal from the noise ¹⁰ ≈4.88*10 ^-4 Hz, the frequency range of the low-frequency coefficient is 0-4.88 x 10 ^-4 Hz, the frequency band occupied by the high-frequency coefficient is 4.88 x 10 ^-4 Hz-0.5 Hz, tidal frequency 9.26 x 10 ^-5 Hz is in the frequency band, when z =11 and z =12 are the same, the tide frequency is also in the frequency band occupied by the decomposed low-frequency coefficient, the tide frequency can be separated from the noise frequency, the tide data can be separated from the noise frequency when the z is more than or equal to 10 and less than or equal to 12 theoretically in the Qingdao semiday tide area, then the tide level sequence based on WGS84 elevation reference can be extracted by setting the high-frequency coefficient to zero and reconstructing the low-frequency coefficient, and the ground tide level sequence is extractedAs shown in fig. 8.

Fig. 3 is a flow chart of the tide level reduction to 1985 national elevation standard based on WGS84 elevation standard, which specifically includes:

r1, acquiring water depth data of a surveying area in a multi-beam sonar strip coverage measurement mode;

r2, performing installation deviation angle correction, sound speed correction and gross error elimination on the water depth data in the step R1, wherein the step can be completed in multi-beam data processing software such as caris;

r3 performs tidal level correction on the corrected water depth data of R2 using the extracted tidal level based on WGS84 elevation reference, and a schematic diagram is shown in fig. 5.

As shown in fig. 6, E represents the position of the survey vessel, and at a certain time E, the elevation anomaly of the grid points P, B, C and D is already solved in step k; as shown in fig. 7, according to the elevation anomalies of the grid points P, B, C and D, a bilinear interpolation algorithm is adopted to solve the elevation anomaly at the ship measuring position;

the quadratic polynomial surface fitting method can also be replaced by the following method: polynomial curve fitting method, cubic spline curve fitting method, three-point fitting method, partition least square plane fitting method, cubic polynomial surface fitting method, hardy multi-face function fitting method, thin plate small-deflection deformation model fitting method, discrete point weighted average fitting method and moving quadratic polynomial fitting method.

It is to be understood that the above description is not intended to limit the present invention, and the present invention is not limited to the above examples, and those skilled in the art may make modifications, alterations, additions or substitutions within the spirit and scope of the present invention.

Claims (8)

1. A method for surveying the altitude of an underway GNSS comprises the following steps:

s1, mounting a measuring device on a ship board, wherein the measuring device comprises a multi-beam sonar, an IMU (inertial measurement unit), a pressure depth sounding sensor and two GNSS (global navigation satellite system) receivers, and the IMU is an inertial measurement unit;

the method comprises the steps that a GNSS receiver independently receives satellite signals, the longitude and latitude and the geodetic elevation of a phase center of the GNSS receiver are obtained, the sonar measures the underwater topography of a measuring area, the longitude and latitude of an underwater measuring point, geodetic height data and the depth are obtained, an IMU obtains a heave value at the IMU and an attitude value of a ship to be measured, the attitude value comprises a rolling value, a pitching value and a yawing value, and a pressure depth measuring sensor measures the distance from the IMU to the water surface in a dynamic draft state when the ship is moving and the distance from the IMU to the water surface when the ship is anchored;

s2, calculating a sea surface geodetic elevation sequence under the WGS84 elevation reference by integrating the data obtained in the S1;

s3, extracting a tide level sequence in the sea surface geodetic elevation sequence under the WGS84 elevation reference by adopting a variational modal decomposition and wavelet transform combined denoising method;

s4, solving the geodetic elevation of the seabed grid points under the WGS84 elevation standard of the measurement area;

s5, obtaining elevation information of 6 or more than 6 GPS leveling combined measuring points near a measuring area, wherein the elevation information simultaneously comprises geodetic elevation of the GPS leveling combined measuring points under WGS84 elevation reference and elevation of the GPS leveling combined measuring points under 1985 national elevation reference;

s6, solving the approximate normal height of the seabed grid points in the measuring area by adopting a quadratic polynomial surface fitting method;

s7, solving the elevation abnormity of the seabed grid points according to the MoroKingsky principle;

and S8, solving out the precise elevation abnormity at the position of the ship to be tested by adopting a bilinear interpolation method according to the longitude and latitude of the phase center of the GNSS receiver in the S1, and reducing the tide level sequence in the S3 to the 1985 national elevation standard to obtain the tide level sequence based on the 1985 national elevation standard.

2. The method for measuring the altitude of the GNSS under navigation according to claim 1, wherein S1 comprises the following substeps:

s1.1, establishing an ideal hull coordinate system by taking an IMU as a reference origin, establishing a three-dimensional model of a measuring device in mapping software, measuring the offset relation between each component on the measuring device through the three-dimensional model, expressing the coordinates of two GNSS receivers under the ideal hull coordinate system, measuring the vertical distance from the reference origin to the GNSS receivers through the three-dimensional model and recording the vertical distance as L1, and obtaining the distance L2 from the IMU to the water surface when a hull is anchored through a pressure depth measuring sensor, namely obtaining the distance L3 from the GNSS receivers to the water surface, wherein the calculation formula is as follows: l3= L1-L2;

s1.2, during sonar measurement, obtaining a roll installation deviation angle and a pitch installation deviation angle of the transducer in a typical area by adopting a 'patch test';

s1.3, correcting the roll value and the pitch value measured by the IMU according to the roll installation deviation angle and the pitch installation deviation angle obtained in the S1.2, and correcting the pitch value measured by the IMU only if the sonar system has a roll compensation function;

s1.4, the accurate geodetic height and longitude and latitude data of the phase centers of the two GNSS receivers are obtained in a GNSS-PPK, GNSS-RTK, network RTK or PPP mode.

3. The method for surveying the altitude of the GNSS, as recited in claim 1, wherein S2 comprises the following substeps:

s2.1, respectively calculating induced heave values H of the two GNSS receivers according to the coordinates of the two GNSS receivers in the S1.1 in the ideal hull coordinate system and the pitch values and the roll values corrected in the S1.3 _i ：H _i ＝xsin P-ysin Rcos P-z(cos Rcos P-1)，H _i Representing induced heave values, x, y and z representing coordinates of the GNSS receiver in an ideal hull coordinate system, and P and R representing pitch values and roll values recorded by the IMU in S1;

s2.2 according to the induced heave value H in S2.1 _i And the heave values at the IMU in S1, respectively calculating the heave values h at the two GNSS receivers _G ：h _G ＝H _i +h ₀ ，h ₀ Represents the heave value at the IMU in S1;

s2.3, the calculation process of the sea surface geodetic elevation sequence is as follows: h = h _GPS -h _G -H _m + delta H, H is the instantaneous ground height of the sea surface, namely the sea surface ground height sequence; h is a total of _GPS The geodetic elevation of the phase center obtained for the GNSS receiver in S1; h _m Is L3; Δ H is the difference between the dynamic draft of the pressure sounding sensor in S1 measuring the motion of the survey vessel and L3 in S1.1.

4. The method for surveying the altitude of the GNSS, as recited in claim 3, wherein S3 comprises the following substeps:

s3.1, successively carrying out K-layer modal decomposition on the sea surface geodetic sequence obtained in the S2 by adopting a variational modal decomposition method, wherein K =10n, n =1,2,3.. Once, n is a positive integer, a penalty coefficient alpha =80000, and a convergence tolerance tol =1e-7;

s3.2, calculating signal energy E obtained by performing K-layer modal decomposition on the sea surface geodetic high sequence each time:

n is the number of sampling points, IMF _K (t) a K-th layer signal sequence for each modal decomposition;

s3.3, solving an energy ratio:

when the energy ratio gamma tends to be stable, the gamma fluctuates at the value of 0, which indicates that K is a proper modal decomposition number;

s3.4, decomposing when the mode decomposition number is K to obtain a first layer mode IMF ₁ Performing variational modal decomposition again, namely performing iterative decomposition for 1 time, referring to S3.2 and S3.3 for the K value judgment criterion of the iterative decomposition for 1 time, and recording the first-layer mode obtained by the iterative decomposition for 1 time as IMF ₁₁ The second layer mode is denoted as IMF ₂₂ IMF of ₁₁ And IMF ₂₂ Reconstructing to obtain a new signal;

s3.5. For pair IMF in S3.4 ₁₁ Performing spectrum analysis to determine IMF ₁₁ The signal spectrum composition of (a);

s3.6. The maximum frequency omega corresponding to the sea surface earth elevation sequence in S2 is as follows:

f _s is the sampling frequency of the GNSS signals;

assuming sea level elevation sequencesFrequency of other elements being f _w The frequency of the tide being f _t The number of layers z using the wavelet number of layers is determined by the following formula:

IMF of S3.4 ₁₁ And IMF ₂₂ Decomposing the z layer by the reconstructed new signal with sym4 wavelet to obtain the z-th layer low-frequency coefficient A _z And a high frequency coefficient D ₁ 、D ₂ ......D _z The frequency band occupied by the low-frequency coefficient is the frequency band of tide, the frequency band occupied by the high-frequency coefficient is the frequency band of other high-frequency elements, the high-frequency coefficient is reset and then reconstructed with the low-frequency coefficient, the forced noise elimination of the sea surface elevation sequence can be realized, and the extracted low-frequency signal is the tide level sequence h based on the WGS84 elevation reference _{tide_84} 。

5. The method for surveying the altitude of the GNSS, as recited in claim 4, wherein S4 comprises the following substeps:

s4.1, carrying out installation deviation angle correction, sound velocity correction and gross error elimination on the water depth data measured by sonar;

s4.2, adopting the h acquired in the S3 to the water depth data processed in the S4.1 _{tide_84} Carrying out tide level correction:

H _i ＝|h _i -h _{tide_84} |，H _i for the elevation of the subsea measurement point, h _i The water depth data after S4.1 treatment is obtained;

s4.3, the sea surface geodetic elevation sequence in the S2 is subjected to grid formation according to 500m × 500m, and geodetic coordinates of grid points under a WGS84 geodetic coordinate system or plane coordinates and geodetic height h are output _Ω Weighting and averaging elevation values of grid points to obtain the average ground height under WGS84 elevation standard, and recording as h _γ 。

6. The method for measuring the altitude of the GNSS under navigation according to claim 5, wherein S6 comprises the following substeps:

s6.1, assuming the existence of the following two times between the elevation abnormity and the plane coordinate on the leveling joint measurement pointPolynomial surface fitting mathematical model:

a ₀ 、a ₁ 、a ₂ 、a ₃ 、a ₄ 、a ₅ for the model to be determined parameter, ζ _i When the number of the leveling joint measuring points is more than 6, the following error equation is formed for the elevation abnormity of the leveling joint measuring points: v = AX- ζ ⁰ A = (X) is obtained according to the least square principle ^T X) ^-1 X ^T ζ ⁰

Solving the undetermined parameter a of the model according to least square primary understanding ₀ 、a ₁ 、a ₂ 、a ₃ 、a ₄ 、a ₅ The value of (d);

s6.2, substituting the coordinates of the seabed grid points in the S4 into the quadratic polynomial surface fitting mathematical model of S6.1 to interpolate approximate elevation abnormity zeta of the seabed grid points ₀ ；

S6.3, measuring the approximate normal height h of the seabed grid points of the area _Γ Comprises the following steps: h is _Γ ＝ζ ₀ +H _i ，ζ ₀ Indicating an approximate elevation anomaly, H _i Is the geodetic elevation of the seabed grid points in S4.

7. The method for measuring the altitude of the GNSS under navigation according to claim 6, wherein S7 comprises the following substeps:

s7.1. Obtaining h obtained in S6 _Γ ；

S7.2, the elevation abnormity zeta of the seabed grid point is formed by a long wave part, namely the approximate elevation abnormity zeta ₀ And short wave part ζ _r Obtaining: ζ = ζ ₀ +ζ _r ；

S7.3, calculating short wave part zeta in elevation abnormity _r ：ζ _r ＝T _r /r，T _r Influence of topographic relief on disturbance positions of the submarine grid points is achieved; zeta _r For terrain correction；

G is universal gravitation constant, G =6.673 × 10 ^-8 C ³ s ^-2 g ^-1 (ii) a ρ is the earth's mean mass density, ρ =2.67gC ^-3 ，γ ₀ ＝[(x-x _p ) ² +(y-y _p ) ² ] ^1/2 (x, y) are the coordinates of the sea-bottom grid point in S4, (x) _p ，y _p ) The coordinate of a certain grid point P in the seabed grid points is obtained;

s7.4, normal gravity r similar to a ground level: r = r ₀ -0.3086h _Γ ，r ₀ To the normal gravity value on the reference ellipsoid, the unit: 10 ^-5 m/s ² ，

Measuring point geographic latitude in a WGS84 geodetic coordinate system in a unit: and (4) degree.

8. The method for surveying the altitude of the GNSS, as recited in claim 7, wherein the S8 comprises the following substeps:

s8.1, according to the longitude and latitude of the phase center of the GNSS receiver in S1, the geodetic elevation of the seabed grid point in S4 and zeta in S7, adopting a bilinear interpolation method to solve zeta which is abnormal in elevation at the position of the ship to be measured ₁ ；

S8.2, calculating the finally interpolated elevation abnormal value f (x, y):

f(x，y ₁ )、f(x，y ₂ ) Showing the elevation anomaly, ζ, interpolated in the x-direction ₁₁ 、ζ ₂₁ 、ζ ₁₂ 、ζ ₂₂ Representing elevation anomaly, x, at four corners of the grid ₁ 、x ₂ 、y ₁ 、y ₂ The coordinates of the four corner points in the X direction and the Y direction are represented, and the coordinates of the ship in the X direction and the Y direction are represented by X and Y;

s8.3, calculating a tide level sequence h based on 1985 national elevation standard _{tide_85} ：h _{tide_85} ＝h _{tide_84} +ζ ₁ 。

Images