CN103605167A - Mallat algorithm-based marine gravity measurement error eliminating method - Google Patents

Abstract

Disclosed in the invention is a Mallat algorithm-based marine gravity measurement error eliminating method. The method comprises the following steps: obtaining latitude, course, and ship speed information outputted by an inertial navigation system as well as a gravity signal measured by a gravity meter in real time; calibrating an initial parameter of the gravity meter; selecting a wavelet function and calculating a scaling function phi (t) and a wavelet function psi (t); calculating coefficients of a high-pass filter and a low-pass filter; carrying out decomposition on the gravity signal by using a Mallat algorithm according to the selected filter coefficients and selecting a decomposition level number; according to a signal to noise ratio of an original signal, calculating a heuristic SURE threshold value and carrying out noise reduction on the gravity signal based on a soft threshold method; reconstructuring the gravity signal after noise reduction, calculating an Etvs correction value by utilizing the information outputted by the inertial navigation system and carrying out filtering processing on the Etvs correction value; and carrying out Etvs correction on the gravity signal after the reconstruction. According to the invention, a defect of insufficient real-time performance during the gravity measurement can be overcome; a main error during the gravity measurement can be eliminated; and the gravity signal precision is improved.

Description

A kind of marine gravimetric survey error cancelling method based on Mallat algorithm

Technical field

The invention belongs to marine gravimetric survey technical field, relate in particular to a kind of marine gravimetric survey error cancelling method based on Mallat algorithm.

Background technology

Last century, the 80-90 age U.S. and USSR (Union of Soviet Socialist Republics) just started to develop the passive navigation backup system of strategic underwater hiding-machine in succession.Initial householder method is based on Graphic Pattern Matching, comprises and the mating of bottom relief map, magnetic chart, but due to needs sonar to measure bottom contour, the disguise that causes submarine topography coupling a little less than; Simultaneously because changes of magnetic field complexity is also difficult to really apply in underwater hiding-machine navigation at present, so gravitational cue and gravity gradient data become the main passive information resources of underwater hiding-machine navigation.Gravity assisting navigation has the advantages such as precision is high, disguised by force, independence is strong, is the desirable positioning means of assisting navigation under water of the submarine navigation devices such as submarine.

The variable signal of gravity is very faint nonstationary random signal, when measuring, is to be all submerged in strong noise background conventionally, how under this strong noise background, to extract gravitational cue value, is the key that improves gravity assisting navigation precision.Traditional Fourier transform processing be stationary signal, for non-stationary signal, need to distinguish various frequency contents, Fourier transform can not meet the requirement that signal is processed.Wavelet analysis is to grow up on the basis of Fourier transform, has made up well the deficiency of Fourier transform.Meanwhile, Mallat algorithm has also guaranteed the real-time that gravitational cue value is obtained.

Summary of the invention

The object of the embodiment of the present invention is to provide a kind of marine gravimetric survey error cancelling method based on Mallat algorithm, what be intended to solve traditional Fourier transform processing is stationary signal, for non-stationary signal, need to distinguish various frequency contents, Fourier transform can not meet the problem of signal processing requirement.

The embodiment of the present invention is achieved in that a kind of marine gravimetric survey error cancelling method based on Mallat algorithm, should comprise the following steps by the marine gravimetric survey error cancelling method based on Mallat algorithm:

Step 1, high-precision inertial navigation system and gravity meter are arranged on same fixed pedestal to the gravitational cue that the information such as latitude, course and the speed of a ship or plane of Real-time Obtaining inertial navigation system output and gravity meter record;

The initial parameter of step 2, gravity meter is demarcated: constant multiplier, zero point drift (mgal/h), the gravity value of gravity reference base station; Before measurement, to demarcate the zero point drift of gravity meter;

Step 3, choose wavelet function, according to the wavelet function of selecting, calculate yardstick function phi (t) and wavelet function ψ (t); If Ψ (t) ∈ is L ²(R), Fourier is changed to Ψ (ω), if met

claim that Ψ (t) is female small echo, stretches and translation to female small echo,

${\mathrm{\Ψ}}_{a,\mathrm{\τ}}\left(t\right)=\frac{1}{\sqrt{a}}\mathrm{\Ψ}\left(\frac{t-\mathrm{\τ}}{a}\right),a>0,\mathrm{\τ}\∈R$

Claim Ψ _{a, τ}(t) be wavelet basis function, wherein a is scale factor or contraction-expansion factor, and τ is called shift factor, and wavelet basis function is acted on to energy signal f (t), obtains continuous wavelet transform:

${\mathrm{WT}}_f(a,\mathrm{\&tau;})=<f\left(t\right),{\mathrm{\&Psi;}}_{a,\mathrm{\&tau;}}\left(t\right)>=\frac{1}{a}{\&Integral;}_{R}f\left(t\right){\mathrm{\&Psi;}}^{*}\left(\frac{t-\mathrm{\&tau;}}{a}\right)\mathrm{dt}$

Wherein, f (t) ∈ L ²(R),, through wavelet transformation, a function of time is projected in yardstick phase plane of two-dimentional time one;

If scaling function is φ (t), corresponding wavelet function is Ψ (t), should meet following two scaling Equations

$\left\{\begin{array}{c}{\mathrm{\&phi;}}_{m,k}\left(t\right)=\underset{l}{\mathrm{\&Sigma;}}{h}_{l-2k}{\mathrm{\&phi;}}_{m+1,l}\left(t\right)\\ {\mathrm{\&Psi;}}_{m,k}=\underset{l}{\mathrm{\&Sigma;}}{g}_{l-2k}{\mathrm{\&phi;}}_{m+1,l}\left(t\right)\end{array}\right.$

Wherein, g (n)=(1) ^1-nh (1-n), definition scaling function is φ (t) ∈ L ²(R), if through stretching on integer translation k and yardstick j, can obtain a yardstick and unique all transformable function sets:

claim that φ (t) is scaling function, h _l-2kand g _l-2kby scaling function φ (t) and wavelet function Ψ (t), determined, be called filter coefficient, φ _{m+1, l}(t) pass through h _l-2kobtain approximation coefficient φ _{m, k}(t), therefore claim h _l-2kfor low-pass filter, φ in like manner _{m+1, l}(t) pass through g _l-2kobtain high frequency details Ψ _{m, k}(t), claim g _l-2kfor Hi-pass filter;

Choose db9 small echo, the support length of db small echo and filter length are all 2N, and vanishing moment is N, the coefficient h of two yardstick difference equations _nmould square have Explicit Expression formula, the support Interval of wavelet function Ψ (t) and scaling function φ (t) is 2N-1, the vanishing moment of Ψ (t) is N; If

wherein

binomial coefficient, h so _ncan be expressed as:

${\left|h_n\right|}^2={\left({\cos }^2\left(\frac{\mathrm{\&omega;}}{2}\right)\right)}^{N}P\left({\sin }^2\left(\frac{\mathrm{\&omega;}}{2}\right)\right)$

Wherein, $m_0\left(\mathrm{\&omega;}\right)=\frac{1}{\sqrt{2}}\underset{k=0}{\overset{2N-1}{\mathrm{\&Sigma;}}}{h}_0{e}^{-\mathrm{ik\&omega;}},$

Step 4, calculating are used for Hi-pass filter and the low-pass filter coefficients of wavelet decomposition; For any f (t) ∈ L ²(R), can be decomposed into the direct sum of several Wavelet Component, so there is the limited decomposed form of following metric space

$V_{m+1}=W_m\&CirclePlus;V_m=W_m\&CirclePlus;W_{m-1}\&CirclePlus;V_{m-1}=\&CenterDot;\&CenterDot;\&CenterDot;=W_m\&CirclePlus;W_{m-1}\&CirclePlus;\&CenterDot;\&CenterDot;\&CenterDot;\&CirclePlus;{\mathrm{<figure-callout\; id="0"\; label="W"\; filenames="HSA0000097532960000042.png"\; state="\{\{state\}\}">W</figure-callout>}}_0\&CirclePlus;{\mathrm{<figure-callout\; id="0"\; label="V"\; filenames="HSA0000097532960000042.png"\; state="\{\{state\}\}">V</figure-callout>}}_0$

Wherein, V _mfor metric space, the metric space being generated by scaling function φ (t), W _mbe called wavelet space, { Ψ _{m, k}(t)=2 ^-m/2Ψ (2 ^-mt-k) } _{k ∈ z, m ∈ z}formed W _morthonormal basis, in decomposable process, each grade of decomposition all can obtain a high frequency details space W _m, be adjacent metric space V _m+1and V _mpoor, due to

be V _m⊥ W _m, therefore, f (t) can use V _mand W _min base jointly launch,

$f\left(t\right)=\underset{k}{\mathrm{\&Sigma;}}{c}_{m,k}{\mathrm{\&phi;}}_{m,k}\left(t\right)+\underset{k}{\mathrm{\&Sigma;}}{d}_{m,k}{\mathrm{\&Psi;}}_{m,k}\left(t\right)$

In formula, the first on the right is the low frequency component of f (t), and second portion is the high fdrequency component of f (t), i.e. the detail section of f (t), and coefficient is wherein

$\left\{\begin{array}{c}{c}_{m,k}=<f\left(t\right),{\mathrm{\&phi;}}_{m,k}\left(t\right)>=\underset{n}{\mathrm{\&Sigma;}}{c}_{m+1,n}\stackrel{\&OverBar;}{h_{m-2k}}\\ d_{m,k}=<f\left(t\right),{\mathrm{\&Psi;}}_{m,k}\left(t\right)>=\underset{n}{\mathrm{\&Sigma;}}{d}_{m+1,n}\stackrel{\&OverBar;}{g_{m-2k}}\end{array}\right.$

From above formula, as known large subspace V _m+1in coefficient c _{m+1, k}time, just can calculate less subspace V _mand W _mcoefficient c _{m, k}and d _{m, k};

Step 5, use Mallat algorithm, decompose gravitational cue according to the filter coefficient having selected, will choose the suitable decomposition number of plies according to different sea status simultaneously; Press formula $f\left(t\right)=\underset{k}{\mathrm{\&Sigma;}}{c}_{m,k}{\mathrm{\&phi;}}_{m,k}\left(t\right)+\underset{k}{\mathrm{\&Sigma;}}{d}_{m,k}{\mathrm{\&Psi;}}_{m,k}\left(t\right)$ Gravitational cue is decomposed;

From above formula, as known large subspace V _m+1in coefficient c _{m+1, k}time, just can calculate less subspace V _mand W _mcoefficient c _{m, k}and d _{m, k};

Step 6, according to the signal to noise ratio (S/N ratio) of original signal, ask for heuristic SURE threshold value and with the method for soft-threshold to gravitational cue noise reduction;

Step 7, to the gravitational cue reconstruct after noise reduction;

Step 8, the information calculating E Tefusi corrected value that utilizes inertial navigation system to export, and E Tefusi corrected value is carried out to filtering processing according to step 3 to step 7;

Step 9, the gravitational cue after reconstruct is carried out to E Tefusi correction.

Further, in step 2: the concrete grammar that the zero point drift of gravity meter is demarcated is:

Suppose on basic point A and B, to have carried out respectively comparison observation when certain marine gravimetric survey starts and finishes, oneself knows that the absolute gravity value of basic point A is g _a, the absolute gravity value that B is ordered is g _b, the difference of the absolute gravity of two basic points is Δ g=g _b-g _a, gravity meter is compared reading and is respectively g on basic point A and B _a' and g _b', difference is Δ g '=g _b'-g _a', the corresponding time of comparison is respectively t _aand t _b, the mistiming is Δ t=t _b-t _a, the zero point drift rate of change of measuring is

C ₀=(Δg-Δg′)/Δt

If measure, start and finish all closures and, in same reference point A, have g _a=g _bso, C ₀can be expressed as

C ₀=(g _A2′-g _A1′)/(t _A2-t _A1)

The difference of supposing the observation date and time of each gravity station that gravity meter duration of work completes and the observation date and time of reference point is followed successively by Δ t ₁, Δ t ₂..., Δ t _n, because the zero point drift modified value of each pendulum point can be calculated as C by linear distribution rule ₀* Δ t _i(i=1,2 ..., n), the zero point drift of gravity meter on i measurement point is modified to

Δg _K=C ₀*Δt _i

The gravity value of each measuring point is g _i=g _i'+Δ g _k=g _i'+C ₀* Δ t _i;

In formula, C ₀unit be mgal/h; Δ t _iunit be h; Δ g _kunit be mgal.

Further, in step 6: choose heuristic SURE threshold value and with the method for soft-threshold to gravitational cue noise reduction, based on Stein, without partial likelihood, estimate that the soft-threshold of SURE estimates it is for given threshold value t, obtains likelihood and estimates, then likelihood function is minimum, obtain required threshold value; Long logarithm threshold value is to be multiplied by a coefficient from obtaining the threshold value t of minimum maximum value variance

the threshold value obtaining;

After trying to achieve threshold value, adopt soft-threshold method that border is occurred to discrete point is retracted to zero, soft-threshold method can effectively be avoided being interrupted, and makes the gravitational cue smoother of rebuilding.

Further, in step 7, to the gravitational cue reconstruct after noise reduction, according to Parseval formula, prove that continuous wavelet exists inverse transformation, inverse transformation formula is

$f\left(t\right)=\frac{1}{C_{\mathrm{\&Psi;}}}{\&Integral;}_0^{+\&infin;}\frac{1}{a^2}\mathrm{da}{\&Integral;}_{-\&infin;}^{+\&infin;}{\mathrm{WT}}_x(a,\mathrm{\&tau;}){\mathrm{\&Psi;}}_{a,\mathrm{\&tau;}}\left(t\right)\mathrm{d\&tau;}$

By selecting wavelet Ψ (t) and to a and τ discretize, obtaining wavelet basis { Ψ _{m, k}(t)=2 ^-m/2Ψ (2 ^-mt-k) } _{k ∈ z, m ∈ z}, the wavelet inverse transformation after discrete can be expressed as

f(t)=∑∑<f(t)，Ψ _m，k(t)>·Ψ _m，k(t)

The reconstruction formula of Mallat algorithm is

$c_{m+1,k}=\underset{n}{\mathrm{\&Sigma;}}{c}_{m,n}{h}_{k-2n}+\underset{n}{\mathrm{\&Sigma;}}{d}_{m,n}{g}_{k-2n}.$

Further, same section of gravitational cue being carried out respectively to three filtering processes.

Further, in step 8, utilize course, position and the velocity information of inertial navigation system output to carry out E Tefusi correction to surveyed gravitational cue, when carrier is when autobiography earth surface moves, the impact that centrifugal force and coriolis force produce the gravity meter being arranged on carrier becomes eotvos effect, the course angle of supposing carrier is ψ, and the speed of a ship or plane is V, and the keel depth of underwater carrier is h; The component velocity of east orientation and north orientation, V _e=Vsin ψ, V _n=Vcos ψ, uses R _eapproximate is the earth radius at L place as latitude, and east orientation component velocity has increased an angular velocity on the basis of earth rotation, and size is

the angular velocity that north orientation ground velocity is corresponding

produce additional centrifugal force

directly act on gravity direction, so, be arranged on gravity that the gravity meter of take on the carrier of speed of a ship or plane V and course ψ experienced as

$g^{\&prime;}=\frac{\mathrm{\&mu;}}{{(R_e-h)}^2}-{(\mathrm{\&omega;}+\frac{V\sin \mathrm{\&psi;}}{(R_e-h)\cos L})}^2(R_e-h){\cos }^{2}L-\frac{V^2{\cos }^2\mathrm{\&psi;}}{R_e-h}$

A/W on earth surface is

$g=\frac{\mathrm{\&mu;}}{{(R_e-h)}^2}-{\mathrm{\&omega;}}^2(R_e-h){\cos }^{2}L$

In the E Tefusi corrected Calculation formula that arrives that is

${\mathrm{\&Delta;g}}_E=g^{\&prime;}-g=2\mathrm{\&omega;}V\sin \mathrm{\&psi;}\cos L+\frac{V^2}{R_e-h}$

The Etvs corrections value of taking h into account differs very little with the Etvs corrections value of ignoring h calculating, the be generally ± 1～± 2mgal of measuring accuracy of current high accuracy gravimeter, and therefore, above formula can be reduced to:

${\mathrm{\&Delta;g}}_E=g^{\&prime;}-g=2\mathrm{\&omega;}V\sin \mathrm{\&psi;}\cos L+\frac{V^2}{R_e}.$

Further, utilize course, latitude and the velocity information of inertial navigation system output all in varying degrees by noise pollution, the E Tefusi corrected value calculating is containing noisy signal, choose same parameter E Tefusi corrected value is also carried out to filtering processing, eliminate the noise that inertial navigation system causes, guarantee that E Tefusi corrected value is corresponding one by one with gravitational cue value in time.

Marine gravimetric survey error cancelling method based on Mallat algorithm provided by the invention, under environment by the ocean complicated, utilize the Mallat algorithm in the present invention to carry out filtering to gravitational cue, and utilize the information of inertial navigation output to carry out E Tefusi correction, compensate to a certain extent the defect of the real-time deficiency in gravimetry, and can eliminate the main error in gravimetry, improve gravitational cue precision.

Accompanying drawing explanation

Fig. 1 is the marine gravimetric survey error cancelling method process flow diagram based on Mallat algorithm that the embodiment of the present invention provides;

Fig. 2 is the gravitational cue value filtering schematic flow sheet that the embodiment of the present invention provides;

Fig. 3 is the Mallat algorithm decomposing schematic representation that the embodiment of the present invention provides;

Fig. 4 is the Mallat algorithm reconstruct schematic diagram that the embodiment of the present invention provides;

Fig. 5 is the actual measurement gravitational cue schematic diagram data that the embodiment of the present invention provides;

Fig. 6 is the filtered gravitational cue schematic diagram data of use Mallat algorithm that the embodiment of the present invention provides;

Fig. 7 is the E Tefusi corrected value schematic diagram that the use inertial navigation data that provide of the embodiment of the present invention are calculated;

Fig. 8 is the filtered E Tefusi corrected value of the use Mallat algorithm schematic diagram that the embodiment of the present invention provides;

Fig. 9 is the gravitational cue schematic diagram data after the filtering of Mallat algorithm and E Tefusi correction that the embodiment of the present invention provides.

Embodiment

In order to make object of the present invention, technical scheme and advantage clearer, below in conjunction with embodiment, the present invention is further elaborated.Should be appreciated that specific embodiment described herein, only in order to explain the present invention, is not intended to limit the present invention.

Below in conjunction with drawings and the specific embodiments, application principle of the present invention is further described.

As shown in Figure 1, the marine gravimetric survey error cancelling method based on Mallat algorithm of the embodiment of the present invention comprises the following steps:

S101: high-precision inertial navigation system and gravity meter are arranged on same fixed pedestal to the gravitational cue that the information such as latitude, course and the speed of a ship or plane of Real-time Obtaining inertial navigation system output and gravity meter record;

S102: the initial parameter of gravity meter is demarcated: constant multiplier, zero point drift (mgal/h), the gravity value of gravity reference base station;

S103: choose suitable wavelet function, according to the wavelet function of selecting, calculate yardstick function phi (t) and wavelet function Ψ (t);

S104: calculate Hi-pass filter and low-pass filter coefficients for wavelet decomposition;

S105: use Mallat algorithm, gravitational cue is decomposed according to the filter coefficient having selected, will choose the suitable decomposition number of plies according to different sea status simultaneously;

S106: according to the signal to noise ratio (S/N ratio) of original signal, ask for heuristic SURE threshold value and with the method for soft-threshold to gravitational cue noise reduction;

S107: to the gravitational cue reconstruct after noise reduction;

S108: utilize the information of inertial navigation system output to calculate E Tefusi corrected value, and E Tefusi corrected value is carried out to filtering processing;

S109: the gravitational cue after reconstruct is carried out to E Tefusi correction.

Concrete steps of the present invention are:

Step 1, first high-precision inertial navigation system and gravity meter are arranged on same fixed pedestal to the gravitational cue that the information such as latitude, course and the speed of a ship or plane of Real-time Obtaining inertial navigation system output and gravity meter record;

Step 2, initialization gravity meter parameter: constant multiplier, zero point drift (mgal/h), the gravity value of gravity reference base station; Variation due to the unstable and environment of the parts of gravimeter own, causes that the reading zero-bit of gravity meter is in continuous variation, therefore before measurement, to demarcate the zero point drift of gravity meter,

Suppose on basic point A and B, to have carried out respectively comparison observation when certain marine gravimetric survey starts and finishes, oneself knows that the absolute gravity value of basic point A is g _a, the absolute gravity value that B is ordered is g _b, the difference of the absolute gravity of two basic points is Δ g=g _b-g _a, gravity meter is compared reading and is respectively g on basic point A and B _a' and g _b', its difference is Δ g '=g _b'-g _a', the corresponding time of comparison is respectively t _aand t _b, its mistiming is Δ t=t _b-t _a, the current zero point drift rate of change of measuring is

C ₀=(Δg-Δg′)/Δt

If measure, start and finish all closures and, in same reference point A, have g _a=g _bso, C ₀can be expressed as

C ₀=(g _A2′-g _A1′)/(t _A2-t _A1)

The difference of supposing the observation date and time of each gravity station that gravity meter duration of work completes and the observation date and time of reference point is followed successively by Δ t ₁, Δ t ₂..., Δ t _n, because the zero point drift modified value of each pendulum point can be calculated as C by linear distribution rule ₀* Δ t _i(i=1,2 ..., n), the zero point drift of gravity meter on i measurement point is modified to

Δg _K=C ₀*Δt _i

The gravity value of each measuring point is g _i=g _i'+Δ g _k=g _i'+C ₀* Δ t _i;

In formula, C ₀unit be mgal/h; Δ t _iunit be h; Δ g _kunit be mgal.

Step 3, choose suitable wavelet function according to the wavelet function of selecting, calculate yardstick function phi (t) and wavelet function Ψ (t), establish Ψ (t) ∈ L ²(R), its Fourier is changed to Ψ (ω), if met

claim that Ψ (t) is female small echo, stretches and translation to female small echo,

${\mathrm{\&Psi;}}_{a,\mathrm{\&tau;}}\left(t\right)=\frac{1}{\sqrt{a}}\mathrm{\&Psi;}\left(\frac{t-\mathrm{\&tau;}}{a}\right),a>0,\mathrm{\&tau;}\&Element;R$

Claim Ψ _{a, τ}(t) be wavelet basis function, wherein a is scale factor or contraction-expansion factor, and τ is called shift factor, and wavelet basis function is acted on to energy signal f (t), obtains continuous wavelet transform:

${\mathrm{WT}}_f(a,\mathrm{\&tau;})=<f\left(t\right),{\mathrm{\&Psi;}}_{a,\mathrm{\&tau;}}\left(t\right)>=\frac{1}{a}{\&Integral;}_{R}f\left(t\right){\mathrm{\&Psi;}}^{*}\left(\frac{t-\mathrm{\&tau;}}{a}\right)\mathrm{dt}$

Wherein, f (t) ∈ L ²(R), through wavelet transformation, a function of time is projected to when two-dimentional in m-yardstick phase plane, be conducive to extract some essential characteristic of signal function;

If scaling function is φ (t), corresponding wavelet function is Ψ (t), and they should meet following two scaling Equations

$\left\{\begin{array}{c}{\mathrm{\&phi;}}_{m,k}\left(t\right)=\underset{l}{\mathrm{\&Sigma;}}{h}_{l-2k}{\mathrm{\&phi;}}_{m+1,l}\left(t\right)\\ {\mathrm{\&Psi;}}_{m,k}=\underset{l}{\mathrm{\&Sigma;}}{g}_{l-2k}{\mathrm{\&phi;}}_{m+1,l}\left(t\right)\end{array}\right.$

Wherein, g (n)=(1) ^1-nh (1-n), definition scaling function is φ (t) ∈ L ²(R), if through stretching on integer translation k and yardstick j, can obtain a yardstick and unique all transformable function sets: claim that φ (t) is scaling function, h _l-2kand g _l-2kby scaling function φ (t) and wavelet function Ψ (t), determined, be conventionally called filter coefficient, φ _{m+1, l}(t) pass through h _l-2kobtain approximation coefficient φ _{m, k}(t), therefore claim h _l-2kfor low-pass filter, φ in like manner _{m+1, l}(t) pass through g _l-2kobtain high frequency details Ψ _{m, k}(t), claim g _l-2kfor Hi-pass filter;

Choose db9 small echo, the support length of db small echo and filter length are all 2N left and right, vanishing moment is N, the wavelet expansion of this sequence is relatively good, can weigh more neatly and increase the boundary problem that support length (in order to put forward high-octane intensity) is brought, although db small echo does not have the expression formula of analytical form, the coefficient h of its pair of yardstick difference equation _nmould square have Explicit Expression formula, the support Interval of wavelet function Ψ (t) and scaling function φ (t) is 2N-1, the vanishing moment of Ψ (t) is N; If

wherein

binomial coefficient, h so _ncan be expressed as:

${\left|h_n\right|}^2={\left({\cos }^2\left(\frac{\mathrm{\&omega;}}{2}\right)\right)}^{N}P\left({\sin }^2\left(\frac{\mathrm{\&omega;}}{2}\right)\right)$

Wherein, $m_0\left(\mathrm{\&omega;}\right)=\frac{1}{\sqrt{2}}\underset{k=0}{\overset{2N-1}{\mathrm{\&Sigma;}}}{h}_0{e}^{-\mathrm{ik\&omega;}};$

Step 4, calculating are used for Hi-pass filter and the low-pass filter coefficients of wavelet decomposition; Calculate the coefficient of dissociation of Hi-pass filter and low-pass filter, for any f (t) ∈ L ²(R), can be decomposed into the direct sum of several Wavelet Component, so there is the limited decomposed form of following metric space

$V_{m+1}=W_m\&CirclePlus;V_m=W_m\&CirclePlus;W_{m-1}\&CirclePlus;V_{m-1}=\&CenterDot;\&CenterDot;\&CenterDot;=W_m\&CirclePlus;W_{m-1}\&CirclePlus;\&CenterDot;\&CenterDot;\&CenterDot;\&CirclePlus;{\mathrm{<figure-callout\; id="0"\; label="W"\; filenames="HSA0000097532960000042.png"\; state="\{\{state\}\}">W</figure-callout>}}_0\&CirclePlus;{\mathrm{<figure-callout\; id="0"\; label="V"\; filenames="HSA0000097532960000042.png"\; state="\{\{state\}\}">V</figure-callout>}}_0$

Wherein, V _mfor metric space,

the metric space being generated by scaling function φ (t), W _mbe called wavelet space, { Ψ _{m, k}(t)=2 ^-m/2Ψ (2 ^-mt-k) } _{k ∈ z, m ∈ z}formed W _morthonormal basis, in decomposable process, each grade of decomposition all can obtain a high frequency details space W _m, it is adjacent metric space V _m+1and V _mpoor, due to

be V _m⊥ W _m, therefore, f (t) can use V _mand W _min base jointly launch,

$f\left(t\right)=\underset{k}{\mathrm{\&Sigma;}}{c}_{m,k}{\mathrm{\&phi;}}_{m,k}\left(t\right)+\underset{k}{\mathrm{\&Sigma;}}{d}_{m,k}{\mathrm{\&Psi;}}_{m,k}\left(t\right)$

In formula, the first on the right is the low frequency component of f (t), and second portion is the high fdrequency component of f (t), i.e. the detail section of f (t), and coefficient is wherein

$\left\{\begin{array}{c}{c}_{m,k}=<f\left(t\right),{\mathrm{\&phi;}}_{m,k}\left(t\right)>=\underset{n}{\mathrm{\&Sigma;}}{c}_{m+1,n}\stackrel{\&OverBar;}{h_{m-2k}}\\ d_{m,k}=<f\left(t\right),{\mathrm{\&Psi;}}_{m,k}\left(t\right)>=\underset{n}{\mathrm{\&Sigma;}}{d}_{m+1,n}\stackrel{\&OverBar;}{g_{m-2k}}\end{array}\right.$

From above formula, as known large subspace V _m+1in coefficient c _{m+1, k}time, just can calculate less subspace V _mand W _mcoefficient c _{m, k}and d _{m, k};

Step 5, as shown in Figure 3, is used Mallat algorithm, and gravitational cue is decomposed according to the filter coefficient having selected, and will choose the suitable decomposition number of plies according to different sea status simultaneously;

Press formula $f\left(t\right)=\underset{k}{\mathrm{\&Sigma;}}{c}_{m,k}{\mathrm{\&phi;}}_{m,k}\left(t\right)+\underset{k}{\mathrm{\&Sigma;}}{d}_{m,k}{\mathrm{\&Psi;}}_{m,k}\left(t\right)$ Gravitational cue is decomposed;

From above formula, as known large subspace V _m+1in coefficient c _{m+1, k}time, just can calculate less subspace V _mand W _mcoefficient c _{m, k}and d _{m, k}, shown in computation process Fig. 3, in Fig. 3, ↓ 2 represent " two extract ", from c _{m+1, k}to d _{m, k}and c _{m, k}, number of samples reduces half, but no matter decomposes how many layers, always the counting and remain unchanged of data, so quadrature discrete wavelet transformation is nonredundant, profit solves the calculated amount of wavelet coefficient in this way lower than the method for direct integration, thereby improves the efficiency that gravitational cue is processed in real time;

Step 6, according to the signal to noise ratio (S/N ratio) of original signal, choose heuristic SURE threshold value and with the method for soft-threshold to gravitational cue noise reduction, choose heuristic SURE threshold value and with the method for soft-threshold to gravitational cue noise reduction, heuristic SURE threshold value is without partial likelihood, to estimate the soft-threshold of (SURE) and the comprehensive form of long logarithm threshold value based on Stein, because the threshold value producing based on SURE threshold value suppresses the DeGrain of noise in the situation that of high signal noise ratio, this method utilizes heuristic function automatically in front kind of threshold value selected, automatically to choose one, wherein, based on Stein, without partial likelihood, estimate the soft-threshold of (SURE) estimates it is for given threshold value t, the likelihood that obtains it is estimated, then likelihood function is minimum, obtain required threshold value, long logarithm threshold value is to be multiplied by a coefficient from obtaining the threshold value t of minimum maximum value variance

the threshold value obtaining,

After trying to achieve threshold value, the method that has two kinds of threshold values on signal, hard threshold method is that the signaling point that the absolute value on a certain decomposition scale is less than to threshold value is thought noise signal, by these signaling point zero setting, it can produce and be interrupted at some point, soft-threshold method is border to be occurred to discrete point is retracted to zero on the basis of hard-threshold, and soft-threshold method can effectively be avoided being interrupted, and makes the gravitational cue smoother of rebuilding;

Step 7, as shown in Figure 4, to the gravitational cue reconstruct after noise reduction, according to Parseval formula, there is inverse transformation in provable continuous wavelet, and its inverse transformation formula is

$f\left(t\right)=\frac{1}{C_{\mathrm{\&Psi;}}}{\&Integral;}_0^{+\&infin;}\frac{1}{a^2}\mathrm{da}{\&Integral;}_{-\&infin;}^{+\&infin;}{\mathrm{WT}}_x(a,\mathrm{\&tau;}){\mathrm{\&Psi;}}_{a,\mathrm{\&tau;}}\left(t\right)\mathrm{d\&tau;}$

By selecting wavelet Ψ (t) and to a and τ discretize, obtaining wavelet basis

wavelet inverse transformation after discrete can be expressed as

f(t)＝∑∑<f(t)，Ψ _m，k(t)>·Ψ _m，k(t)

The reconstruction formula of Mallat algorithm is

$c_{m+1,k}=\underset{n}{\mathrm{\&Sigma;}}{c}_{m,n}{h}_{k-2n}+\underset{n}{\mathrm{\&Sigma;}}{d}_{m,n}{g}_{k-2n}$

Restructuring procedure Fig. 2

Simultaneously, for relatively and the operation efficiency of the operation efficiency algorithm of other wavelet algorithms, same section of gravitational cue carried out respectively to three filtering to be processed, contrast the working time to common wavelet algorithm and Wavelet Packet Algorithm, known according to table 1, Mallat algorithm operation efficiency is better than other wavelet algorithms;

Table 1Mallat algorithm and Wavelet Packet Algorithm comparison working time (unit: s)

Step 8, the information calculating E Tefusi corrected value that utilizes inertial navigation system to export, and E Tefusi corrected value is carried out to filtering processing by step 3 to step 7; Utilize course, position and the velocity information of inertial navigation system output to carry out E Tefusi correction to surveyed gravitational cue, when carrier is when autobiography earth surface moves, the impact that centrifugal force and coriolis force produce the gravity meter being arranged on carrier becomes eotvos effect, the course angle of supposing carrier is ψ, the speed of a ship or plane is V, and the keel depth of underwater carrier is h; The component velocity of east orientation and north orientation, V _e=Vsin ψ, V _n=Vcos ψ, uses R _eapproximate is the earth radius at L place as latitude, and east orientation component velocity has increased an angular velocity on the basis of earth rotation, and size is

the angular velocity that north orientation ground velocity is corresponding

produce additional centrifugal force

directly act on gravity direction, so, be arranged on gravity that the gravity meter of take on the carrier of speed of a ship or plane V and course ψ experienced as

$g^{\&prime;}=\frac{\mathrm{\&mu;}}{{(R_e-h)}^2}-{(\mathrm{\&omega;}+\frac{V\sin \mathrm{\&psi;}}{(R_e-h)\cos L})}^2(R_e-h){\cos }^{2}L-\frac{V^2{\cos }^2\mathrm{\&psi;}}{R_e-h}$

A/W on earth surface is

$g=\frac{\mathrm{\&mu;}}{{(R_e-h)}^2}-{\mathrm{\&omega;}}^2(R_e-h){\cos }^{2}L$

In the E Tefusi corrected Calculation formula that arrives that is

${\mathrm{\&Delta;g}}_E=g^{\&prime;}-g=2\mathrm{\&omega;}V\sin \mathrm{\&psi;}\cos L+\frac{V^2}{R_e-h}$

The Etvs corrections value of taking h into account differs very little with the Etvs corrections value of ignoring h calculating, the be generally ± 1～± 2mgal of measuring accuracy of current high accuracy gravimeter, and therefore, above formula can be reduced to:

${\mathrm{\&Delta;g}}_E=g^{\&prime;}-g=2\mathrm{\&omega;}V\sin \mathrm{\&psi;}\cos L+\frac{V^2}{R_e}$

Utilize the course of inertial navigation system output, latitude and velocity information are all in varying degrees by noise pollution, the E Tefusi corrected value calculating is also containing noisy signal, if directly gravitational cue value is proofreaied and correct, precision will reduce, secondly, in processing the filtering of gravitational cue, output valve is effect after the filtering of data in a period of time, the measured value gathering is not the gravity value on a certain independent point, but the reflection of mean gravity value in a period of time, but corresponding each value of putting the separately during output of inertial navigation, therefore above-mentioned E Tefusi corrected value is all also single point in time, if remove adjusted weight force signal by these values, on time, can not accomplish complete correspondence, error also just increases thereupon, choose same parameter E Tefusi corrected value is also carried out to filtering processing, not only can eliminate the noise that inertial navigation system causes, can also guarantee that E Tefusi corrected value is corresponding one by one with gravitational cue value in time,

Step 9, the gravitational cue after reconstruct is carried out to E Tefusi correction.

Principle of work of the present invention:

The present invention is for general gravitational cue acquisition methods, because the cycle of low-pass filter is longer, and E Tefusi to proofread and correct be all that the data that provide by GPS are carried out post-processed, this causes the gravitational cue cannot Real-time Obtaining.Use Mallat algorithm to carry out filtering processing, and in conjunction with the information of inertial navigation system output, gravitational cue is carried out to E Tefusi correction, improved the efficiency of gravitational cue data processing, make precision have certain guarantee simultaneously, not only can be for gravity assisting navigation and seafari, and can also carry out the correction of Gravity Models on this basis.Therefore, the research of the marine gravimetric survey error cancelling method based on Mallat algorithm has great importance.

In conjunction with following experiment, excellent beneficial effect of the present invention is described further:

The original signal that Fig. 5 measures for testing gravity meter on certain river, gravity meter and optical fibre gyro inertial navigation system are arranged on same pedestal, usage data receives the signal that software receives gravity meter and inertial navigation system output simultaneously, guarantee temporal synchronous, in experiment, effect for checking E Tefusi proofreaies and correct, has the turning that is repeatedly greater than 30 °, accompanying drawing 5～9th, gravitational cue result in the process of moving.

From Fig. 5, find out, the noise of raw data is very strong, and signal to noise ratio (S/N ratio) is very low, useful gravitational cue is submerged completely, and the peak value of interference noise is up to 350mgal/h, and choosing of visible wave filter is the key of removing noise, Fig. 6 is through the filtered gravitational cue of Mallat algorithm, through contrast, find, Mallat algorithm has effectively suppressed noise, recovers useful gravitational cue.

The foregoing is only preferred embodiment of the present invention, not in order to limit the present invention, all any modifications of doing within the spirit and principles in the present invention, be equal to and replace and improvement etc., within all should being included in protection scope of the present invention.

Claims (7)

1. the marine gravimetric survey error cancelling method based on Mallat algorithm, is characterized in that, should comprise the following steps by the marine gravimetric survey error cancelling method based on Mallat algorithm:

Step 1, high-precision inertial navigation system and gravity meter are arranged on same fixed pedestal to the gravitational cue that latitude, course and the speed of a ship or plane information of Real-time Obtaining inertial navigation system output and gravity meter record;

The initial parameter of step 2, gravity meter is demarcated: constant multiplier, zero point drift, the gravity value of gravity reference base station; Before measurement, to demarcate the zero point drift of gravity meter;

Step 3, choose wavelet function, according to the wavelet function of selecting, calculate yardstick function phi (t) and wavelet function Ψ (t); If Ψ (t) ∈ is L ²(R), Fourier is changed to Ψ (ω), if met

claim that Ψ (t) is female small echo, stretches and translation to female small echo,

${\mathrm{\&Psi;}}_{a,\mathrm{\&tau;}}\left(t\right)=\frac{1}{\sqrt{a}}\mathrm{\&Psi;}\left(\frac{t-\mathrm{\&tau;}}{a}\right),a>0,\mathrm{\&tau;}\&Element;R$

Claim Ψ _{a, τ}(t) be wavelet basis function, wherein a is scale factor or contraction-expansion factor, and τ is called shift factor, and wavelet basis function is acted on to energy signal f (t), obtains continuous wavelet transform:

${\mathrm{WT}}_f(a,\mathrm{\&tau;})=<f\left(t\right),{\mathrm{\&Psi;}}_{a,\mathrm{\&tau;}}\left(t\right)>=\frac{1}{a}{\&Integral;}_{R}f\left(t\right){\mathrm{\&Psi;}}^{*}\left(\frac{t-\mathrm{\&tau;}}{a}\right)\mathrm{dt}$

Wherein, f (t) ∈ L ²(R), through wavelet transformation, a function of time is projected to when two-dimentional in m-yardstick phase plane;

If scaling function is φ (t), corresponding wavelet function is Ψ (t), should meet following two scaling Equations

$\left\{\begin{array}{c}{\mathrm{\&phi;}}_{m,k}\left(t\right)=\underset{l}{\mathrm{\&Sigma;}}{h}_{l-2k}{\mathrm{\&phi;}}_{m+1,l}\left(t\right)\\ {\mathrm{\&Psi;}}_{m,k}=\underset{l}{\mathrm{\&Sigma;}}{g}_{l-2k}{\mathrm{\&phi;}}_{m+1,l}\left(t\right)\end{array}\right.$

Wherein, g (n)=(1) ^1-nh (1-n), definition scaling function is φ (t) ∈ L ²(R), if through stretching on integer translation k and yardstick j, can obtain a yardstick and unique all transformable function sets:

claim that φ (t) is scaling function, h _l-2kand g _l-2kby scaling function φ (t) and wavelet function Ψ (t), determined, be called filter coefficient, φ _{m+1, l}(t) pass through h _l-2kobtain approximation coefficient φ _{m, k}(t), therefore claim h _l-2kfor low-pass filter, φ in like manner _{m+1, l}(t) pass through g _l-2kobtain high frequency details Ψ _{m, k}(t), claim g _l-2kfor Hi-pass filter;

Choose db9 small echo, the support length of db small echo and filter length are all 2N, and vanishing moment is N, the coefficient h of two yardstick difference equations _nmould square have Explicit Expression formula, the support Interval of wavelet function Ψ (t) and scaling function φ (t) is 2N-1, the vanishing moment of Ψ (t) is N; If

wherein binomial coefficient, h so _ncan be expressed as:

${\left|h_n\right|}^2={\left({\cos }^2\left(\frac{\mathrm{\&omega;}}{2}\right)\right)}^{N}P\left({\sin }^2\left(\frac{\mathrm{\&omega;}}{2}\right)\right)$

Wherein, $m_0\left(\mathrm{\&omega;}\right)=\frac{1}{\sqrt{2}}\underset{k=0}{\overset{2N-1}{\mathrm{\&Sigma;}}}{h}_0{e}^{-\mathrm{ik\&omega;}},$

Step 4, calculating are used for Hi-pass filter and the low-pass filter coefficients of wavelet decomposition; For any f (t) ∈ L ²(R), can be decomposed into the direct sum of several Wavelet Component, so there is the limited decomposed form of following metric space

$V_{m+1}=W_m\&CirclePlus;V_m=W_m\&CirclePlus;W_{m-1}\&CirclePlus;V_{m-1}=\&CenterDot;\&CenterDot;\&CenterDot;=W_m\&CirclePlus;W_{m-1}\&CirclePlus;\&CenterDot;\&CenterDot;\&CenterDot;\&CirclePlus;W_0\&CirclePlus;V_0$

Wherein, V _mfor metric space

the metric space being generated by scaling function φ (t), W _mbe called wavelet space, { Ψ _{m, k}(t)=2 ^-m/2Ψ (2 ^-mt-k) } _{k ∈ z, m ∈ z}formed W _morthonormal basis, in decomposable process, each grade of decomposition all can obtain a high frequency details space W _m, be adjacent metric space V _m+1and V _mpoor, due to

be V _m⊥ W _m, therefore, f (t) can use V _mand W _min base jointly launch,

$f\left(t\right)=\underset{k}{\mathrm{\&Sigma;}}{c}_{m,k}{\mathrm{\&phi;}}_{m,k}\left(t\right)+\underset{k}{\mathrm{\&Sigma;}}{d}_{m,k}{\mathrm{\&Psi;}}_{m,k}\left(t\right)$

In formula, the first on the right is the low frequency component of f (t), and second portion is the high fdrequency component of f (t), i.e. the detail section of f (t), and coefficient is wherein

$\left\{\begin{array}{c}{c}_{m,k}=<f\left(t\right),{\mathrm{\&phi;}}_{m,k}\left(t\right)>=\underset{n}{\mathrm{\&Sigma;}}{c}_{m+1,n}\stackrel{\&OverBar;}{h_{m-2k}}\\ d_{m,k}=<f\left(t\right),{\mathrm{\&Psi;}}_{m,k}\left(t\right)>=\underset{n}{\mathrm{\&Sigma;}}{d}_{m+1,n}\stackrel{\&OverBar;}{g_{m-2k}}\end{array}\right.$

From above formula, as known large subspace V _m+1in coefficient c _{m+1, k}time, just can calculate less subspace V _mand W _mcoefficient c _{m, k}and d _{m, k};

Step 5, use Mallat algorithm, decompose gravitational cue according to the filter coefficient having selected, will choose the decomposition number of plies according to different sea status simultaneously; Press formula

gravitational cue is decomposed;

From above formula, as known large subspace V _m+1in coefficient c _{m+1, k}time, just can calculate less subspace V _mand W _mcoefficient c _{m, k}and d _{m, k};

Step 6, according to the signal to noise ratio (S/N ratio) of original signal, ask for heuristic SURE threshold value and with the method for soft-threshold to gravitational cue noise reduction;

Step 7, to the gravitational cue reconstruct after noise reduction;

Step 8, the information calculating E Tefusi corrected value that utilizes inertial navigation system to export, and E Tefusi corrected value is carried out to filtering processing according to step 3 to step 7;

Step 9, the gravitational cue after reconstruct is carried out to E Tefusi correction.

2. the marine gravimetric survey error cancelling method based on Mallat algorithm as claimed in claim 1, is characterized in that, in step 2: the concrete grammar that the zero point drift of gravity meter is demarcated is:

Suppose on basic point A and B, to have carried out respectively comparison observation when certain marine gravimetric survey starts and finishes, oneself knows that the absolute gravity value of basic point A is g _a, the absolute gravity value that B is ordered is g _b, the difference of the absolute gravity of two basic points is Δ g=g _b-g _a, gravity meter is compared reading and is respectively g on basic point A and B _a' and g _b', difference is Δ g '=g _b'-g _a', the corresponding time of comparison is respectively t _aand t _b, the mistiming is Δ t=t _b-t _a, the zero point drift rate of change of measuring is

C ₀＝(Δg-Δg′)/Δt

If measure, start and finish all closures and, in same reference point A, have g _a=g _bso, C ₀can be expressed as

C ₀=(g _A2′-g _A1′)/(t _A2-t _A1)

The difference of supposing the observation date and time of each gravity station that gravity meter duration of work completes and the observation date and time of reference point is followed successively by Δ t ₁, Δ t ₂..., Δ t _n, because the zero point drift modified value of each pendulum point can be calculated as C by linear distribution rule ₀* Δ t _i(i=1,2 ..., n), the zero point drift of gravity meter on i measurement point is modified to

Δg _K=C ₀*Δt _i

The gravity value of each measuring point is g _i=g _i'+Δ g _k=g _i'+C ₀* Δ t _i;

In formula, C ₀unit be mgal/h; Δ t _iunit be h; Δ g _kunit be mgal.

3. the marine gravimetric survey error cancelling method based on Mallat algorithm as claimed in claim 1, it is characterized in that, in step 6: choose heuristic SURE threshold value and with the method for soft-threshold to gravitational cue noise reduction, based on Stein, without partial likelihood, estimate the soft-threshold of SURE estimates it is for given threshold value t, obtaining likelihood estimates, then likelihood function is minimum, obtain required threshold value; Long logarithm threshold value is to be multiplied by a coefficient from obtaining the threshold value t of minimum maximum value variance

the threshold value obtaining;

After trying to achieve threshold value, adopt soft-threshold method that border is occurred to discrete point is retracted to zero, soft-threshold method can effectively be avoided being interrupted, and makes the gravitational cue smoother of rebuilding.

4. the marine gravimetric survey error cancelling method based on Mallat algorithm as claimed in claim 1, is characterized in that, in step 7, to the gravitational cue reconstruct after noise reduction, according to Parseval formula, proves that continuous wavelet exists inverse transformation, and inverse transformation formula is

$f\left(t\right)=\frac{1}{C_{\mathrm{\&Psi;}}}{\&Integral;}_0^{+\&infin;}\frac{1}{a^2}\mathrm{da}{\&Integral;}_{-\&infin;}^{+\&infin;}{\mathrm{WT}}_x(a,\mathrm{\&tau;}){\mathrm{\&Psi;}}_{a,\mathrm{\&tau;}}\left(t\right)\mathrm{d\&tau;}$

By selecting wavelet Ψ (t) and to a and τ discretize, obtaining wavelet basis { Ψ _{m, k}(t)=2 ^-m/2Ψ (2 ^-mt-k) } _{k ∈ z, m ∈ z}, the wavelet inverse transformation after discrete can be expressed as

f(t)=∑∑<f(t)，Ψ _m，k(t)>·Ψ _m，k(t)

The reconstruction formula of Mallat algorithm is

$c_{m+1,k}=\underset{n}{\mathrm{\&Sigma;}}{c}_{m,n}{h}_{k-2n}+\underset{n}{\mathrm{\&Sigma;}}{d}_{m,n}{g}_{k-2n}.$

5. the marine gravimetric survey error cancelling method based on Mallat algorithm as claimed in claim 4, is characterized in that, same section of gravitational cue is carried out respectively to three filtering and process.

6. the marine gravimetric survey error cancelling method based on Mallat algorithm as claimed in claim 1, it is characterized in that, in step 8, utilize course, position and the velocity information of inertial navigation system output to carry out E Tefusi correction to surveyed gravitational cue, when carrier is when autobiography earth surface moves, the impact that centrifugal force and coriolis force produce the gravity meter being arranged on carrier becomes eotvos effect, and the course angle of supposing carrier is ψ, the speed of a ship or plane is V, and the keel depth of underwater carrier is h; The component velocity of east orientation and north orientation, V _e=Vsin ψ, V _n=Vcos ψ, uses R _eapproximate is the earth radius at L place as latitude, and east orientation component velocity has increased an angular velocity on the basis of earth rotation, and size is the angular velocity that north orientation ground velocity is corresponding

produce additional centrifugal force

directly act on gravity direction, so, be arranged on gravity that the gravity meter of take on the carrier of speed of a ship or plane V and course ψ experienced as

$g^{\&prime;}=\frac{\mathrm{\&mu;}}{{(R_e-h)}^2}-{(\mathrm{\&omega;}+\frac{V\sin \mathrm{\&psi;}}{(R_e-h)\cos L})}^2(R_e-h){\cos }^{2}L-\frac{V^2{\cos }^2\mathrm{\&psi;}}{R_e-h}$

A/W on earth surface is

$g=\frac{\mathrm{\&mu;}}{{(R_e-h)}^2}-{\mathrm{\&omega;}}^2(R_e-h){\cos }^{2}L$

In the E Tefusi corrected Calculation formula that arrives that is

${\mathrm{\&Delta;g}}_E=g^{\&prime;}-g=2\mathrm{\&omega;}V\sin \mathrm{\&psi;}\cos L+\frac{V^2}{R_e-h}$

The Etvs corrections value of taking h into account differs very little with the Etvs corrections value of ignoring h calculating, the be generally ± 1～± 2mgal of measuring accuracy of current high accuracy gravimeter, and therefore, above formula can be reduced to:

${\mathrm{\&Delta;g}}_E=g^{\&prime;}-g=2\mathrm{\&omega;}V\sin \mathrm{\&psi;}\cos L+\frac{V^2}{R_e}.$

7. the marine gravimetric survey error cancelling method based on Mallat algorithm as claimed in claim 6, it is characterized in that, utilize course, latitude and the velocity information of inertial navigation system output all in varying degrees by noise pollution, the E Tefusi corrected value calculating is containing noisy signal, choose same parameter E Tefusi corrected value is also carried out to filtering processing, eliminate the noise that inertial navigation system causes, guarantee that E Tefusi corrected value is corresponding one by one with gravitational cue value in time.

Images