CN109856690A - Aeromagnetic tensor data noise restraint method and system based on mixing norm fitting - Google Patents
Aeromagnetic tensor data noise restraint method and system based on mixing norm fitting Download PDFInfo
- Publication number
- CN109856690A CN109856690A CN201910150563.1A CN201910150563A CN109856690A CN 109856690 A CN109856690 A CN 109856690A CN 201910150563 A CN201910150563 A CN 201910150563A CN 109856690 A CN109856690 A CN 109856690A
- Authority
- CN
- China
- Prior art keywords
- data
- compensation
- signal
- change
- time rate
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Abstract
The present invention discloses a kind of aeromagnetic tensor data noise restraint method and system based on mixing norm fitting.Method includes: to establish magnetic gradient aggregation of data compensation model；According to the magnetic gradient comprehensive compensation model, the time rate of change signal of observation data is obtained；According to the time rate of change signal, the compensation equation based on change rate signal is constructed；According to the compensation equation based on change rate signal, data fit object function is established；The fit object function is solved, compensating parameter is obtained；According to the compensating parameter, compensated magnetic gradient data are determined；Compensation effect is assessed according to the compensated magnetic gradient data.It can be improved the signal-to-noise ratio of final thermal compensation signal using method or apparatus of the invention.
Description
Technical field
The present invention relates to geophysics airborne magnetic prospecting technical fields, more particularly to one kind based on mixing norm fitting
Aeromagnetic tensor data noise restraint method and system.
Background technique
Airborne magnetic prospecting is a kind of geophysical exploration technology hand that is important and being widely used in mineral products and environment survey
Section.In recent years, the magnetic gradient measurements system of based superconductive quantum inteferometer (SQUID) is greatly developed, so that directly measuring high-precision
The magnetic gradient data of degree are possibly realized, and the full tensor magnetic gradiometer of superconduction is because of its high information quantity, highfield sensitivity, small in size etc.
Plurality of advantages, it is considered to be the developing direction of third generation boat magnetic detection.Currently, China is built to be equipped with the full tensor based on SQUID
Magnetic gradient measurements platform prototype, different from traditional aeromagnetic gradient survey flying platform measured using airborne mode,
Magnetic gradient measurements instrument and relevant device (GPS, inertial navigation system (INS), reading reading circuit etc.) are loaded into one by the prototype
In gondola, and flight measurement is carried out in such a way that helicopter pulls, which will substantially reduce helicopter and interfere magnetic measurement
Influence, or even pull hawser it is sufficiently long under the conditions of helicopter influence can be ignored, but there are still big in gondola
The equipment that can generate magnetic disturbance is measured, therefore effectively handles the measurement accuracy for improving flying platform by magnetic compensation, to realize height
Resolution ratio magnetic gradient measurements have great importance.
Since magnetic gradient Tensor measuring system uses the highly sensitive characteristic of SQUID sensor, in low latitude, prospecting flies
Influence of the geomagnetic gradient for sensor cannot be ignored in row, and the compensation system for utilizing traditional high-altitude test flight to obtain
It is several, magnetic gradient signal in low-latitude flying can not be compensated for the interference effect of sensor, therefore directly use traditional compensation side
Method (i.e. by high-altitude flight data obtain penalty coefficient directly low latitude observation flight data are compensated in the way of) can not obtain
Obtain the magnetic gradient processing data that precision is met the requirements.Therefore it is directed to this problem, needs to change original compensating operation mode, directly
It connects and low-latitude flying survey data is compensated.But the major issue for needing to face during this is, induction generates interference source
Magnetic gradient signal simultaneously be also likely to be for prospecting purpose useful signal, then how suppress interference under the premise of retain
These useful signals improve the signal-to-noise ratio of final reconstruction signal, become an important problem of magnetic gradient Data processing.
Summary of the invention
The object of the present invention is to provide a kind of aeromagnetic tensor data noise restraint method based on mixing norm fitting and it is
System, can be improved the signal-to-noise ratio of final thermal compensation signal.
To achieve the above object, the present invention provides following schemes:
A kind of aeromagnetic tensor data noise restraint method based on mixing norm fitting, comprising:
Establish magnetic gradient aggregation of data compensation model；
According to the magnetic gradient comprehensive compensation model, the time rate of change signal of observation data is obtained；
According to the time rate of change signal, the compensation equation based on change rate signal is constructed；
According to the compensation equation based on change rate signal, data fit object function is established；
The fit object function is solved, compensating parameter is obtained；
According to the compensating parameter, compensated magnetic gradient data are determined；
Compensation effect is assessed according to the compensated magnetic gradient data.
Optionally, described according to the magnetic gradient comprehensive compensation model, obtain the time rate of change signal of observation data, tool
Body includes:
According to the magnetic gradient comprehensive compensation model, building compensation equation Ax=b；
According to the compensation equation, time rate of change the signal DA, Db of observation data are obtained；
Wherein, A is the observation data directly acquired, and b is the observation data for compensating gradient value, and x is in practical calculate
The compensating parameter used, DA are the time rate of change signal of the observation data directly acquired, and Db is the sight for compensating gradient value
The time rate of change signal of measured data.
Optionally, described that the compensation equation based on change rate signal is constructed according to the time rate of change signal, it is specific to wrap
It includes:
According to the time rate of change signal, the compensation equation DAx=Db based on change rate signal is constructed；
Wherein, x is compensating parameter used in practical calculate, and DA is the time rate of change letter of the observation data directly acquired
Number, Db is the time rate of change signal for compensating the observation data of gradient value.
Optionally, described that the fit object function is solved, compensating parameter is obtained, is specifically included:
Based on Lp norm data fitting method, initial compensating parameter and initial compensation data are obtained；
Calculate the single order time rate of change of the initial compensation data；
The provincial characteristics that useful signal is analyzed according to the single order time rate of change obtains analysis result；
The value of Lp norm order p and q is determined according to the analysis result；
According to the value of the p and q, the data fit object is solved by weighting Gauss-Newton iterative algorithm again
Function obtains compensating parameter.
Optionally, described that compensation effect is assessed according to the compensated magnetic gradient data, it specifically includes:
According to the compensated magnetic gradient data, compensation effect is assessed using normalizing parameter, the standard
Changing parameter includes root mean square, standard deviation and improvement ratio.
A kind of aeromagnetic tensor data noise suppression system based on mixing norm fitting, comprising:
Compensation model establishes module, for establishing magnetic gradient aggregation of data compensation model；
Time rate of change signal acquisition module, for obtaining observation data according to the magnetic gradient comprehensive compensation model
Time rate of change signal；
Equation determining module is compensated, for constructing the compensation based on change rate signal according to the time rate of change signal
Equation；
Objective function determining module, for establishing data fitting mesh according to the compensation equation based on change rate signal
Scalar functions；
It solves module and obtains compensating parameter for solving to the fit object function；
Magnetic gradient data determining module, for determining compensated magnetic gradient data according to the compensating parameter；
Evaluation module, for being assessed according to the compensated magnetic gradient data compensation effect.
Optionally, the time rate of change signal acquisition module, specifically includes:
First compensation equation determination unit, for according to the magnetic gradient comprehensive compensation model, building compensation equation Ax=
b；
Time rate of change signal acquiring unit, for obtaining the time rate of change of observation data according to the compensation equation
Signal DA, Db；
Wherein, A is the observation data directly acquired, and b is the observation data for compensating gradient value, and x is in practical calculate
The compensating parameter used, DA are the time rate of change signal of the observation data directly acquired, and Db is the sight for compensating gradient value
The time rate of change signal of measured data.
Optionally, the compensation equation determining module, specifically includes:
Second compensation equation determination unit, for constructing based on change rate signal according to the time rate of change signal
Compensate equation DAx=Db；
Wherein, x is compensating parameter used in practical calculate, and DA is the time rate of change letter of the observation data directly acquired
Number, Db is the time rate of change signal for compensating the observation data of gradient value.
Optionally, the solution module, specifically includes:
Based on Lp norm data fitting method, initial compensating parameter and initial compensation data are obtained；
Calculate the single order time rate of change of the initial compensation data；
The provincial characteristics that useful signal is analyzed according to the single order time rate of change obtains analysis result；
The value of Lp norm order p and q is determined according to the analysis result；
According to the value of the p and q, adaptive Lp norm data are solved by weighting Gauss-Newton iterative algorithm again
The objective function of fitting.
Optionally, the evaluation module, specifically includes:
Assessment unit, for being carried out to compensation effect using normalizing parameter according to the compensated magnetic gradient data
Assessment, the normalizing parameter include root mean square, standard deviation and improvement ratio.
The specific embodiment provided according to the present invention, the invention discloses following technical effects: the present invention provides a kind of base
In the aeromagnetic tensor data noise restraint method of mixing norm fitting, the letter of magnetic gradient signal after final compensation can effectively improve
It makes an uproar and compares, backoff algorithm is enabled to facilitate the compensation data for being directly suitable for highly sensitive magnetic gradient sensor.
Detailed description of the invention
It in order to more clearly explain the embodiment of the invention or the technical proposal in the existing technology, below will be to institute in embodiment
Attached drawing to be used is needed to be briefly described, it should be apparent that, the accompanying drawings in the following description is only some implementations of the invention
Example, for those of ordinary skill in the art, without any creative labor, can also be according to these attached drawings
Obtain other attached drawings.
Fig. 1 is that the present invention is based on the aeromagnetic tensor data noise restraint method flow charts of mixing norm fitting；
Fig. 2 is magnetic susceptibility Distribution Value figure of the present invention；
Fig. 3 is that magnetic gradient sensor G1 channel measurement signal of the present invention and distinct methods compensation data values compare；
Fig. 4 is that magnetic gradient sensor G2 channel measurement signal of the present invention and distinct methods compensation data values compare；
Fig. 5 is that magnetic gradient sensor G3 channel measurement signal of the present invention and distinct methods compensation data values compare；
Fig. 6 is that magnetic gradient sensor G4 channel measurement signal of the present invention and distinct methods compensation data values compare；
Fig. 7 is that magnetic gradient sensor G5 channel measurement signal of the present invention and distinct methods compensation data values compare；
Fig. 8 is that magnetic gradient sensor G6 channel measurement signal of the present invention and distinct methods compensation data values compare；
Fig. 9 is that the present invention is based on the aeromagnetic tensor data noise suppression system construction drawings of mixing norm fitting；
Algorithm shows the compensation of the gradometer measured value in the practical channel G1 and filtering processing result to Figure 10 through the invention
It is intended to, wherein (a) indicates the preceding comparing result with data after compensation of compensation；(b) indicate compensate processing after, filtering before with
The comparing result of filtered data；
After Figure 11 is compensated and is filtered to the gradometer measured value in the channel G1, data after measurement data, compensation,
The Fourier Frequency and Amplitude spectrum analysis comparison diagram of filtered data.
Specific embodiment
Following will be combined with the drawings in the embodiments of the present invention, and technical solution in the embodiment of the present invention carries out clear, complete
Site preparation description, it is clear that described embodiments are only a part of the embodiments of the present invention, instead of all the embodiments.It is based on
Embodiment in the present invention, it is obtained by those of ordinary skill in the art without making creative efforts every other
Embodiment shall fall within the protection scope of the present invention.
The object of the present invention is to provide a kind of aeromagnetic tensor data noise restraint method based on mixing norm fitting and it is
System, can be improved the signal-to-noise ratio of final thermal compensation signal.
In order to make the foregoing objectives, features and advantages of the present invention clearer and more comprehensible, with reference to the accompanying drawing and specific real
Applying mode, the present invention is described in further detail.
Fig. 1 is that the present invention is based on the aeromagnetic tensor data noise restraint method flow charts of mixing norm fitting.Such as Fig. 1 institute
Show, a kind of aeromagnetic tensor data noise restraint method based on mixing norm fitting, comprising:
Step 101: establishing magnetic gradient aggregation of data compensation model；
According to generality it is assumed that the signal of gradiometer measurement is by actual signal, low frequency interfering noise and high frequency electromagnetic noise
Composition, wherein low frequency interfering noise is mainly derived from external magnetic source interference (remanent magnetism, the sense magnetic, eddy current magnetism of measuring table), with
And not calibrated remaining sensing system deviation and stray magnetic fields etc..Therefore, when consider measuring table interference source mainly from
It inside gondola, and is when being rigidly connected with Magnetic Sensor, the measurement gradient value of Magnetic Sensor can carry out table by such as drag
Show:
Wherein, l is Measurement channel number；For the magnetic gradient value of measurement；For the true magnetic ladder at channel measurement point
Angle value；HeFor local background magnetic field, can be obtained by Global Geomagnetic Field Models IGRF-12 model or regional geomagnetic field model are approximate
It arrives；dHe=(dHx,dHy,dHz)=(He(t+1/fs)-He(t)) fs is background magnetic field time rate of change signal, and wherein fs is
Data sampling frequency；BrFor the three axis magnetometer measured value after compensation correction；For not calibrated remaining system deviation,It is outer
The remanent magnetism of portion's interference source interferes, and the two is invariable normal due to corresponding to single gradiometer measured value within a certain period of time
Value, can integrate as a deviation factor ol；For the linear operator of external interference source sense magnetic disturbance, by (2p+1) × 3 coefficient
Composition；For the linear operator of external interference source vortex interference, it is made of (2p+1) × 3 coefficient；h0For stray magnetic fields interference
Linear operator, be made of (2q+1) × 3 coefficient；N is high frequency Gauss white noise.By to above-mentioned acquisition signal model into
Row arranges, and is available for the computation model of magnetic gradient compensation data, as follows:
Wherein, A is by all kinds of observation data and relevant calculation amountThe N of composition
× M rank matrix, N are gradiometer sampling number；For reality
Compensating parameter used in the calculating of border, total M=12p+6q+10.Prolong since the model considers the system between each sensor
Slow characteristic (inertial navigation system, gradiometer and magnetometer), thus it is dry to the outside that different sensors measurement data is calculated is related to
The sense magnetic, vortex interference field and stray magnetic fields interference for disturbing source carry out linear filtering approximate processing；When do not consider system delay spy
Property when, can use p=q=0, then model can simplify as with M=10 penalty coefficient x=(ol,ax,ay,az,bx,by,bz,cx,
cz,cz)TComputation model, can indicate are as follows:
Step 102: according to the magnetic gradient comprehensive compensation model, obtaining the time rate of change signal of observation data；
According to the magnetic gradient comprehensive compensation model, building compensation equation Ax=b；
According to the compensation equation, time rate of change the signal DA, Db of observation data are obtained；
Wherein, A is the observation data directly acquired, and b is the observation data for compensating gradient value, and x is in practical calculate
The compensating parameter used, DA are the time rate of change signal of the observation data directly acquired, and Db is the sight for compensating gradient value
The time rate of change signal of measured data.
Step 103: according to the time rate of change signal, constructing the compensation equation based on change rate signal；It specifically includes:
According to the time rate of change signal, the compensation equation DAx=Db based on change rate signal is constructed；
Wherein, x is compensating parameter used in practical calculate, and DA is the time rate of change letter of the observation data directly acquired
Number, Db is the time rate of change signal for compensating the observation data of gradient value.
Step 104: according to the compensation equation based on change rate signal, establishing data fit object function；
Fig. 2 is magnetic susceptibility Distribution Value figure of the present invention.According to theoretical model calculate it is found that if target signal source be it is multiple not
Magnetic source with morphological feature forms, then a series of (earth coordinates-measuring table coordinate system-single ladder after being converted by coordinates
Spend flowmeter sensor and measure parallel direction), magnetic anomaly has form similar to magnetic gradient abnormal signal to it under earth coordinates
Property, only there is small magnitude relevant to flight attitude variation to fluctuate in the biggish middle section of magnetic anomaly amplitude.Therefore magnetic anomaly letter
It number can be divided into two kinds of regions substantially, the first kind is the flat section of relative smooth in non-magnetic source region or remote magnetic source region, with
And second class be magnetic source area or magnetic source borderline region magnetic fluctuation area.In order to keep the noise with higher of signal after compensating
Than needing to reduce the invalid fluctuation of signal, the especially low-frequency fluctuation in signal as far as possible in first kind region, because of high frequency
Fluctuation can then need to protect as far as possible by low-pass filtering treatment without influencing last compensation effect in the second class region
Hold small signal concussion relevant to flight attitude variation in the fluctuation amplitude and form and region of signal.
In general sense, the mode of traditional L2 norm data fitting, as it is assumed that the mistake between calculated value and observed quantity
Difference meets Gaussian distribution, therefore can be very good to be fitted interference signal when useful signal amplitude is very low, and
So that final thermal compensation signal is intended to small magnitude and uniformly shakes or constant value；But when calculated value and the direct residual volume of observed quantity
(Gaussian distribution is unsatisfactory for) when with the variation of obvious amplitude, since this method can carry out indifference fitting to useful signal,
Data are in the apparent invalid signals concussion of useful signal areas adjacent appearance after finally making compensation, thus reduce thermal compensation signal
Signal-to-noise ratio.And the fitting of Lp norm (0≤p≤1) data not only can effectively keep the step evolution of signal, moreover it is possible to reduce wherein
The invalid concussion of reconstruction signal boundary.In addition to this, the derivative information of signal can further display the borderline properties of signal,
Thus according to above-mentioned analysis, in conjunction with the time-derivative information of Lp norm and signal, following objective function is constructed:
Wherein, p, q are the order of Lp norm,OrThe p that each region uses, q value can be adaptively chosen according to the characteristic of signal in practical calculating；
α is that fitting coefficient of balance can generally take for balancing the weight of two parts fitting data
Due to the objective function non-differentiability in the case of p, q < 2, therefore for ease of calculation usually by its above-mentioned objective function
Discrete form is approximately a kind of objective function form of weighted L2 norm, i.e., are as follows:
Wherein, Φ (x) is objective function, W1(x),W2It (x) is norm weighting matrix, diagonal element w1ii,w2ii, specifically
Definition mode is shown in following formula:
And there is εp,εqIt is the sufficiently small constant of a value, can is 1e-8. with value for when calculating in practical calculating
Between change come gradient operator matrix be expressed as form:
The above problem can be iterated solution by weight Weighted Gauss-newton Gauss-Newton method.
Step 105: the fit object function being solved, compensating parameter is obtained；
Step 105 specifically includes:
Based on Lp norm data fitting method, initial compensating parameter and initial compensation data are obtained；
Calculate the single order time rate of change of the initial compensation data；
The provincial characteristics that useful signal is analyzed according to the single order time rate of change obtains analysis result；
The value of Lp norm order p and q is determined according to the analysis result；
According to the value of the p and q, the data fit object is solved by weighting Gauss-Newton iterative algorithm again
Function obtains compensating parameter.
The specific steps of step 105 are as follows:
(1) step 1:
Since the fit object function adopts measurement data and measurement data change rate in different data areas
It is fitted norm order with different Lp, to improve the signal-to-noise ratio of final thermal compensation signal, it is therefore desirable to an initial offset data
It is inputted as algorithm, to judge the initial possibility region where magnetic anomaly source.It is fitted in this method using general L2 norm
As a result initial compensation data are provided, ridge regression algorithm can be used or the method for solving for being based on principal component analysis (PCA) in method for solving,
Specific numerical procedure is as follows:
(1) the data fit object function based on L2 norm are as follows:
(2) pass through the available corresponding compensating parameter x of method of ridge regressionβ, and offset data is calculatedBase
This calculation method are as follows:
Wherein, β is ridge regression parameter, and I is M rank unit matrix.
(3) by the method for solving based on PCA, calculated when sytem matrix linear dependence is more serious it is more stable,
It is specific to calculate step are as follows:
(i) coefficient matrices A is carried out to mean value standardization by column is
Wherein, aiFor the column vector of matrix A,For vector aiMean value, N is number of sampling points.
(ii) rightCarry out Eigenvalues Decomposition C=V Σ VT, obtain eigenvectors matrix V=(v1,v2,K,vM),
And corresponding eigenvalue λ1> λ2> K > λM, (Σ=diag (λ1,λ2,K,λM))
(iii) selected characteristic value is greater than λminComponent m (≤M) it is a, constitute principal component transform penalty coefficient matrix P=
AVm,Vm=(v1,v2,K,vm)；Wherein, P is the principal component matrix of A, Vm=(v1,v2,K,vm) it is by preceding m feature vector components makes
The submatrix of composition, viFor the ith feature vector in character pair vector matrix V.
(iv) data G after solving principal component transform compensating parameter y by least square method and compensating accordinglyrl0:
(2) step 2:
Data after the compensation obtained according to initial L2 norm data fitting relevant calculationCalculating judges demarcation signal
In non-signal region, smooth signal area and Non-smooth surface signal area, circular it is as follows:
(1) first to signal after compensationAnd its rate of change valueIt is normalized, obtains normalization data
di,bi:
(2) normalized signal d, the threshold θ of b are set1,θ2, and according to following rule judgement and value data and its variation
The order p, q of rate fitting norm:
The parameter will be used in following data the Fitting Calculation.
(3) step 3:
According to order p, the q of the fitting norm that the above method determines, base is iteratively solved by weighting Gauss-Newton again
In the objective function of Lp mixing norm data fitting, specific export process is as follows.
In order to acquire the minimum of objective function, need to meet target function gradient value equal to 0, due to former objective function
In weighted term W1(x),W2(x) related to compensating parameter x, therefore can be unfolded to obtain the close of former problem by single order Taylor
Like objective function, and may further obtain iterating to calculate formula accordingly:
According to condition existing for above-mentioned objective function minimal point:
Following iterative calculation formula can be further derived by:
Wherein, the coefficient matrix in the linear iteraction formula and right-hand vector can be embodied as:
In formulaFor norm weighting matrix, diagonal element w1ii,w2ii, mode is specifically defined as following formula institute
Show:
And there is εp,εqIt is the sufficiently small constant of a value, can is 1e-8 with value in practical calculating.
Therefore according to above-mentioned iterative calculation formula, data can be by calculating as follows after compensating parameter x and accordingly compensation
Method is calculated:
(i) initial value x is provided0, iteration maximum step number maxstep, and stagnate error e rror；
(ii) doi=1 ..., maxstep；
(ii-1) weighting coefficient matrix is updated
(ii-2) Gauss-Hessian matrix H is constructednAnd gradient vector gn:
(ii-3) by linear iterative equation group derivation algorithm (such as PCG/GMRES/BiCG) or direct method, (LU is decomposed
Deng) solve system of equation HnΔ x=-gn, obtain updating variable Δ x, and update calculated value xn+1/2=xn+Δx；
(ii-4) according to constraint condition, constraint section projection x is carried out to more new variables is calculatedn+1=Π |Q(xn+1/2)；
(ii-5) examine whether iteration is stagnated, if met | | xn+1-xn||/||xn| | < error then stagnate by iteration, otherwise
(ii-1) is gone to continue to calculate.
(iii) data after compensating are calculated by compensating parameter x
It should be noted that the data volume and data class for constructing magnetic compensation model equation coefficient matrices A in calculating
Type difference (related to the modeling pattern of interferer models), the coefficient matrices A of formation can have serious pathosis (i.e. conditional number
It is very big), in order to improve the reusability of the stability and compensating parameter that solve parameter, it can use the PCA in above-mentioned steps 1 (3)
The principal component matrix P of analysis method extraction coefficient matrix A, and construct data fit object function in the same way and solve
Data after to compensating parameter y and compensationIt is as follows that objectives function and offset data solve step:
Step 106: according to the compensating parameter, determining compensated magnetic gradient data；
Step 107: compensation effect being assessed according to the compensated magnetic gradient data.
Generally data after compensation evaluate compensation result using following normalizing parameter: root mean square RMS, standard deviation
σc, improve compare IR.Being defined as follows for parameter value is shown:
Wherein, dcFor the data after compensation deals, d0For original acquisition data, N is same type of data amount check.
In the case where carrying out theoretical model test, due to known accurate signal value dt, therefore can be by believing after compensation
Residual r number with accurate signalc=dt-dcThe validity for carrying out quantization algorithm, be defined as follows three evaluation parameters: residual is equal
Root RMS, residual standard deviation sigmac, Signal to Noise Ratio (SNR) and Y-PSNR PSNR, shown in being defined as follows of parameter value:
Wherein, residual root mean square RMS, residual standard deviation sigmacFor measuring parametric error, it is worth the smaller the better；Signal-to-noise ratio
SNR and Y-PSNR PSNR is used to measure the reconstruction accuracy of useful signal, and value is the bigger the better.
Fig. 9 is that the present invention is based on the aeromagnetic tensor data noise suppression system construction drawings of mixing norm fitting.Such as Fig. 9 institute
Show, a kind of aeromagnetic tensor data noise suppression system based on mixing norm fitting, comprising:
Compensation model establishes module 201, for establishing magnetic gradient aggregation of data compensation model；
Time rate of change signal acquisition module 202, for obtaining observation data according to the magnetic gradient comprehensive compensation model
Time rate of change signal；
Equation determining module 203 is compensated, for constructing the benefit based on change rate signal according to the time rate of change signal
Repay equation；
Objective function determining module 204, for establishing data fitting according to the compensation equation based on change rate signal
Objective function；
It solves module 205 and obtains compensating parameter for solving to the fit object function；
Magnetic gradient data determining module 206, for determining compensated magnetic gradient data according to the compensating parameter；
Evaluation module 207, for being assessed according to the compensated magnetic gradient data compensation effect.
The time rate of change signal acquisition module 202, specifically includes:
First compensation equation determination unit, for according to the magnetic gradient comprehensive compensation model, building compensation equation Ax=
b；
Time rate of change signal acquiring unit, for obtaining the time rate of change of observation data according to the compensation equation
Signal DA, Db；
Wherein, A is the observation data directly acquired, and b is the observation data for compensating gradient value, and x is in practical calculate
The compensating parameter used, DA are the time rate of change signal of the observation data directly acquired, and Db is the sight for compensating gradient value
The time rate of change signal of measured data.
The compensation equation determining module 203, specifically includes:
Second compensation equation determination unit, for constructing based on change rate signal according to the time rate of change signal
Compensate equation DAx=Db；
Wherein, x is compensating parameter used in practical calculate, and DA is the time rate of change letter of the observation data directly acquired
Number, Db is the time rate of change signal for compensating the observation data of gradient value.
The solution module 205, specifically includes:
Based on Lp norm data fitting method, initial compensating parameter and initial compensation data are obtained；
Calculate the single order time rate of change of the initial compensation data；
The provincial characteristics that useful signal is analyzed according to the single order time rate of change obtains analysis result；
The value of Lp norm order p and q is determined according to the analysis result.
According to the value of the p and q, adaptive Lp norm data are solved by weighting Gauss-Newton iterative algorithm again
The objective function of fitting.
The evaluation module 207, specifically includes:
Assessment unit, for being carried out to compensation effect using normalizing parameter according to the compensated magnetic gradient data
Assessment, the normalizing parameter include root mean square, standard deviation and improvement ratio.
Specific embodiment 1:
The effect of institute's magnetic gradient compensation data method of the present invention is carried out below by theoretical model data and measured data
Verifying and explanation.SQUID magnetic gradient measurements instrument involved in the present invention shares 6 magnetic gradient sensor output channels: G1-G6.
Wherein, in measured data verifying, the present embodiment says magnetic gradient compensation data effect using the data of G1 channel measurement
It is bright.
The effect of magnetic gradient compensation data method proposed by the invention is verified and said by emulating data first
It is bright.It assuming that the magnetic susceptibility value of basement rock is very small, is 1e-8 (SI) to be embedded in two magnetic sources that magnetic susceptibility value is different, form is different
Body, left side plate body are heterogeneous source, and magnetic susceptibility value is followed successively by 1e-3,2e-3,3e-3 (SI) from top to bottom, and right side is hung down
The magnetic susceptibility value of straight tube-like body is 5e-3 (SI), as shown in Figure 3.It is true by inertial navigation system and GPS under certain observed altitude
The posture position information of measurement, magnetic anomaly signal coordinates is transformed under flight measurement coordinate system, obtains approximate the earth background magnetic
Field signal F0, magnetic anomaly field three-component Fv and magnetic gradient component signal Gi, set magnetic anomaly constant value as Fv, then corresponding local magnetic field
Data can be expressed as Fa=F0+Fv.Assuming that the interference source of measuring table is equivalent to unit dipole, the surplus of hypothesis is utilized
Magnetic, sense magnetic and vortex parameter and degree of unbalancedness coefficient, then can be constructed by following model construction comprising air maneuver noise
Actual measurement gradient value:
Wherein, T be earth coordinates arrive sensor measurement plane coordinate conversion matrix, used in for construct benefit
The data for repaying model coefficient matrix include: F0, dF0.Pass through three kinds of data fit object function method for solving: ridge regression, PCA divides
Analysis, the method for the present invention obtains compensated magnetic gradient data such as Fig. 5, and shown in 6, wherein Fig. 3 is magnetic gradient sensor G1 of the present invention
Channel measurement signal and distinct methods compensation data values compare.Different curves respectively indicate what distinct methods were calculated in Fig. 3
Offset data, Fig. 4 are that magnetic gradient sensor G2 channel measurement signal of the present invention and distinct methods compensation data values compare；Fig. 5 is
Magnetic gradient sensor G3 channel measurement signal of the present invention and distinct methods compensation data values compare；Fig. 6 is magnetic gradient of the present invention biography
Sensor G4 channel measurement signal and distinct methods compensation data values compare；Fig. 7 is magnetic gradient sensor G5 channel measurement of the present invention
Signal and distinct methods compensation data values compare；Fig. 8 is magnetic gradient sensor G6 channel measurement signal of the present invention and distinct methods
Compensation data values compare；All kinds of evaluation indexes of data are as shown in table 1- table 3 after corresponding compensation.Pass through series of computation result
It has been shown that, new method can effectively keep the amplitude and minutia of useful signal compared to conventional method, and evaluation parameter is each
Entry value is intended to be substantially better than other two methods.
The processing data evaluation parametric statistics of magnetic gradient compensation method of the table 1 based on ridge regression data fitting solution scheme
The processing data evaluation parametric statistics of magnetic gradient compensation method of the table 2 based on PCA analysis data fitting solution scheme
The processing data evaluation parametric statistics of magnetic gradient compensation method of the table 3 based on new method data fitting solution scheme
Then, it is verified and is illustrated by effect of the measured data to institute's magnetic gradient compensation data method of the present invention.It adopts
Sample data height is 1km, therefore does not include effective magnetic gradient signal.Algorithm is terraced to the magnetic in the practical channel G1 through the invention by Figure 10
The compensation and filtering processing result schematic diagram of degree measurement magnitude, wherein (a) indicates the comparison knot before compensating with data after compensation
Fruit；(b) it indicates after compensating processing, the preceding comparing result with filtered data of filtering.Gradometer of the Figure 11 to the channel G1
After measured value is compensated and is filtered, data, the Fourier Frequency and Amplitude of filtered data spectrum after measurement data, compensation
Analyze comparison diagram.By Figure 10 and Figure 11 it is found that the magnetic compensation method can effectively suppress the low-frequency disturbance in measurement data,
And offset data improves compensation precision by can further eliminate high-frequency noise therein after filtering processing.Correspondingly, compensation
Quality evaluation index corresponding with filtering processing result is provided by table 4, in conjunction with Fig. 9~11 and table 4 it is found that through overcompensation and filtering
Treated, and magnetic gradient data root mean square is less than 20pT, improves than IR up to 2.0e3, therefore is based on method energy proposed by the present invention
Good compensation result is enough obtained, can reach the requirement of instrumental resolution.
Processing data evaluation parametric statistics of the table 4 based on comprehensive magnetic gradient compensation method
So far, attached drawing is had been combined basic skills and embodiment of the invention is described in detail.Foundation is retouched above
It states, those skilled in the art should be to the magnetic gradient number of the invention based on the single acquisition channel signal of SQUID aeromagnetic gradiometer
There is clear understanding according to compensation method.It should be noted that in attached drawing or specification text, the realization that is not painted or describes
Mode is form known to a person of ordinary skill in the art in technical field, is not described in detail.
Each embodiment in this specification is described in a progressive manner, the highlights of each of the examples are with other
The difference of embodiment, the same or similar parts in each embodiment may refer to each other.For system disclosed in embodiment
For, since it is corresponded to the methods disclosed in the examples, so being described relatively simple, related place is said referring to method part
It is bright.
Used herein a specific example illustrates the principle and implementation of the invention, and above embodiments are said
It is bright to be merely used to help understand method and its core concept of the invention；At the same time, for those skilled in the art, foundation
Thought of the invention, there will be changes in the specific implementation manner and application range.In conclusion the content of the present specification is not
It is interpreted as limitation of the present invention.
Claims (10)
1. a kind of aeromagnetic tensor data noise restraint method based on mixing norm fitting characterized by comprising
Establish magnetic gradient aggregation of data compensation model；
According to the magnetic gradient comprehensive compensation model, the time rate of change signal of observation data is obtained；
According to the time rate of change signal, the compensation equation based on change rate signal is constructed；
According to the compensation equation based on change rate signal, data fit object function is established；
The fit object function is solved, compensating parameter is obtained；
According to the compensating parameter, compensated magnetic gradient data are determined；
Compensation effect is assessed according to the compensated magnetic gradient data.
2. the aeromagnetic tensor data noise restraint method according to claim 1 based on mixing norm fitting, feature exist
In, it is described according to the magnetic gradient comprehensive compensation model, the time rate of change signal of observation data is obtained, is specifically included:
According to the magnetic gradient comprehensive compensation model, building compensation equation Ax=b；
According to the compensation equation, time rate of change the signal DA, Db of observation data are obtained；
Wherein, A is the observation data directly acquired, and b is the observation data for compensating gradient value, and x is to use in practical calculate
Compensating parameter, DA is the time rate of change signal of the observation data directly acquired, and Db is the observation number for compensating gradient value
According to time rate of change signal.
3. the aeromagnetic tensor data noise restraint method according to claim 1 based on mixing norm fitting, feature exist
In, it is described according to the time rate of change signal, the compensation equation based on change rate signal is constructed, is specifically included:
According to the time rate of change signal, the compensation equation DAx=Db based on change rate signal is constructed；
Wherein, x is compensating parameter used in practical calculate, and DA is the time rate of change signal of the observation data directly acquired,
Db is the time rate of change signal for compensating the observation data of gradient value.
4. the aeromagnetic tensor data noise restraint method according to claim 1 based on mixing norm fitting, feature exist
In, it is described that the fit object function is solved, compensating parameter is obtained, is specifically included:
Based on Lp norm data fitting method, initial compensating parameter and initial compensation data are obtained；
Calculate the single order time rate of change of the initial compensation data；
The provincial characteristics that useful signal is analyzed according to the single order time rate of change obtains analysis result；
The value of fitting Lp norm order p and q is determined according to the analysis result；
According to the value of the p and q, the data fit object letter is solved by weighting Gauss-Newton iterative algorithm again
Number, obtains compensating parameter.
5. the boat magnetic tensor data noise restraint method according to claim 4 based on extension compensation model, which is characterized in that institute
It states and compensation effect is assessed according to the compensated magnetic gradient data, specifically include:
According to the compensated magnetic gradient data, compensation effect is assessed using normalizing parameter, the standardization ginseng
Number includes root mean square, standard deviation and improvement ratio.
6. a kind of aeromagnetic tensor data noise suppression system based on mixing norm fitting characterized by comprising
Compensation model establishes module, for establishing magnetic gradient aggregation of data compensation model；
Time rate of change signal acquisition module, for obtaining the time of observation data according to the magnetic gradient comprehensive compensation model
Change rate signal；
Equation determining module is compensated, for constructing the compensation equation based on change rate signal according to the time rate of change signal；
Objective function determining module, for establishing data fit object letter according to the compensation equation based on change rate signal
Number；
It solves module and obtains compensating parameter for solving to the fit object function；
Magnetic gradient data determining module, for determining compensated magnetic gradient data according to the compensating parameter；
Evaluation module, for being assessed according to the compensated magnetic gradient data compensation effect.
7. the aeromagnetic tensor data noise suppression system according to claim 6 based on mixing norm fitting, feature exist
In the time rate of change signal acquisition module specifically includes:
First compensation equation determination unit, for according to the magnetic gradient comprehensive compensation model, building compensation equation Ax=b；
Time rate of change signal acquiring unit, for obtaining the time rate of change signal of observation data according to the compensation equation
DA,Db；
Wherein, A is the observation data directly acquired, and b is the observation data for compensating gradient value, and x is to use in practical calculate
Compensating parameter, DA is the time rate of change signal of the observation data directly acquired, and Db is the observation number for compensating gradient value
According to time rate of change signal.
8. the aeromagnetic tensor data noise suppression system according to claim 6 based on mixing norm fitting, feature exist
In the compensation equation determining module specifically includes:
Second compensation equation determination unit, for constructing the compensation based on change rate signal according to the time rate of change signal
Equation DAx=Db；
Wherein, x is compensating parameter used in practical calculate, and DA is the time rate of change signal of the observation data directly acquired,
Db is the time rate of change signal for compensating the observation data of gradient value.
9. the aeromagnetic tensor data noise suppression system according to claim 6 based on mixing norm fitting, feature exist
In the solution module specifically includes:
Based on Lp norm data fitting method, initial compensating parameter and initial compensation data are obtained；
Calculate the single order time rate of change of the initial compensation data；
The provincial characteristics that useful signal is analyzed according to the single order time rate of change obtains analysis result；
The value of Lp norm order p and q is determined according to the analysis result；
According to the value of the p and q, adaptive Lp norm data fitting is solved by weighting Gauss-Newton iterative algorithm again
Objective function.
10. the boat magnetic tensor data noise suppression system according to claim 6 based on extension compensation model, which is characterized in that
The evaluation module, specifically includes:
Assessment unit, for being assessed compensation effect using normalizing parameter according to the compensated magnetic gradient data,
The normalizing parameter includes root mean square, standard deviation and improvement ratio.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910150563.1A CN109856690B (en) | 2019-02-28 | 2019-02-28 | Aeromagnetic gradient tensor data noise suppression method and system based on mixed norm fitting |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910150563.1A CN109856690B (en) | 2019-02-28 | 2019-02-28 | Aeromagnetic gradient tensor data noise suppression method and system based on mixed norm fitting |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109856690A true CN109856690A (en) | 2019-06-07 |
| CN109856690B CN109856690B (en) | 2020-04-28 |
Family
ID=66899267
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910150563.1A Expired - Fee Related CN109856690B (en) | 2019-02-28 | 2019-02-28 | Aeromagnetic gradient tensor data noise suppression method and system based on mixed norm fitting |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109856690B (en) |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111079285A (en) * | 2019-12-16 | 2020-04-28 | 山东大学 | Full-tensor magnetic gradient data compensation optimization method and system |
| CN112666619A (en) * | 2020-12-17 | 2021-04-16 | 中国自然资源航空物探遥感中心 | Method and device for obtaining accurate aviation magnetic measurement data based on standard deviation method |
| CN113420423A (en) * | 2021-06-01 | 2021-09-21 | 郑州步始智能科技有限公司 | Method for forward modeling full tensor from magnetic field vertical component based on regularization method in frequency domain |
| CN113534292A (en) * | 2021-07-09 | 2021-10-22 | 吉林大学 | Small signal model aeromagnetic compensation method based on forgetting factor RLS |
| CN115391726A (en) * | 2022-07-26 | 2022-11-25 | 睿励科学仪器（上海）有限公司 | Fitting optimization method for ellipsometry system and related device |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106959471A (en) * | 2017-04-21 | 2017-07-18 | 中国科学院电子学研究所 | Boat magnetic compensation method based on the non-linear boat total field gradient compensation model of magnetic |
| CN107291659A (en) * | 2017-05-16 | 2017-10-24 | 哈尔滨工程大学 | The recurrence cosine transform method of the step upward continuation plane modulus gradient of plane GEOMAGNETIC FIELD one |
| CN108227024A (en) * | 2017-12-04 | 2018-06-29 | 中国科学院地质与地球物理研究所 | A kind of method and system using full tensor magnetic gradient data inversion underground magnetic susceptibility |
-
2019
- 2019-02-28 CN CN201910150563.1A patent/CN109856690B/en not_active Expired - Fee Related
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106959471A (en) * | 2017-04-21 | 2017-07-18 | 中国科学院电子学研究所 | Boat magnetic compensation method based on the non-linear boat total field gradient compensation model of magnetic |
| CN107291659A (en) * | 2017-05-16 | 2017-10-24 | 哈尔滨工程大学 | The recurrence cosine transform method of the step upward continuation plane modulus gradient of plane GEOMAGNETIC FIELD one |
| CN108227024A (en) * | 2017-12-04 | 2018-06-29 | 中国科学院地质与地球物理研究所 | A kind of method and system using full tensor magnetic gradient data inversion underground magnetic susceptibility |
Non-Patent Citations (1)
| Title |
|---|
| SHUANGXI JI ET AL.: "Compensation of SQUID gradient data in magnetic tensor gradiometer system for airborne exploration", 《THE INTERNATIONAL SYMPOSIUM ON DEEP EARTH EXPLORATION AND PRACTICES》 * |
Cited By (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111079285A (en) * | 2019-12-16 | 2020-04-28 | 山东大学 | Full-tensor magnetic gradient data compensation optimization method and system |
| CN111079285B (en) * | 2019-12-16 | 2022-02-08 | 山东大学 | Full-tensor magnetic gradient data compensation optimization method and system |
| CN112666619A (en) * | 2020-12-17 | 2021-04-16 | 中国自然资源航空物探遥感中心 | Method and device for obtaining accurate aviation magnetic measurement data based on standard deviation method |
| CN112666619B (en) * | 2020-12-17 | 2021-08-24 | 中国自然资源航空物探遥感中心 | Method and device for obtaining accurate aviation magnetic measurement data based on standard deviation method |
| CN113420423A (en) * | 2021-06-01 | 2021-09-21 | 郑州步始智能科技有限公司 | Method for forward modeling full tensor from magnetic field vertical component based on regularization method in frequency domain |
| CN113534292A (en) * | 2021-07-09 | 2021-10-22 | 吉林大学 | Small signal model aeromagnetic compensation method based on forgetting factor RLS |
| CN113534292B (en) * | 2021-07-09 | 2022-12-27 | 吉林大学 | Small signal model aeromagnetic compensation method based on forgetting factor RLS |
| CN115391726A (en) * | 2022-07-26 | 2022-11-25 | 睿励科学仪器（上海）有限公司 | Fitting optimization method for ellipsometry system and related device |
| CN115391726B (en) * | 2022-07-26 | 2023-10-27 | 睿励科学仪器（上海）有限公司 | Fitting optimization method and related device for ellipsometry system |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109856690B (en) | 2020-04-28 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109856690A (en) | Aeromagnetic tensor data noise restraint method and system based on mixing norm fitting | |
| CN109814163B (en) | Method and system for suppressing noise of aeromagnetic tensor data based on extended compensation model | |
| CN110095738B (en) | Parameter design method for magnetic gradiometer | |
| CN111433634B (en) | Magnetic compensation method based on aeromagnetic compensation error model | |
| CN104931028B (en) | A kind of three axle magneto-electronic compass error compensation methods based on deep learning | |
| Han et al. | A modified Tolles–Lawson model robust to the errors of the three-axis strapdown magnetometer | |
| CN107389049B (en) | Magnetic compass error real-time compensation method based on quasi-Kalman factor | |
| Hardwick | Important design considerations for inboard airborne magnetic gradiometers | |
| CN109856689B (en) | Noise suppression processing method and system for superconducting aeromagnetic gradient tensor data | |
| AU2008231589B2 (en) | Terrain correction systems | |
| Dou et al. | An aeromagnetic compensation coefficient-estimating method robust to geomagnetic gradient | |
| CN107132587B (en) | The full tensor magnetic gradient measurements system mounting error calibration method of aviation superconduction and device | |
| CN113074721B (en) | Geomagnetic fingerprint construction method based on magnetic moment method | |
| CN110187393B (en) | Aeromagnetic compensation method based on generalized regression neural network | |
| CN110274586A (en) | Magnetic air compensation method comprising the compensation of more spectrum atom magnetometer deflection errors | |
| CN110018429B (en) | Method and system for eliminating magnetic field interference caused by vibration of magnetic detection platform | |
| CN112858961B (en) | Compensation method for multisource magnetic interference of aviation platform | |
| CN103562753A (en) | Gravity gradiometer survey techniques | |
| CN108520112A (en) | A kind of aircraft interference field compensation method based on Tikhonov regularization | |
| CN108344409A (en) | Method for improving satellite attitude determination precision | |
| Wang et al. | Compensation for mobile carrier magnetic interference in a SQUID-based full-tensor magnetic gradiometer using the flower pollination algorithm | |
| CN108254796B (en) | A kind of optimization method of scalar magnetometer array baseline | |
| Mu et al. | Small Signal Magnetic Compensation Method for UAV-Borne Vector Magnetometer System | |
| Yin | Compensation for aircraft effects of magnetic tensor gradiometer with nonlinear least square method | |
| CN114325848A (en) | Aviation magnetic interference compensation method and device for adaptive correction of compensation coefficient |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20200428 |