CN109856690B - Aeromagnetic gradient tensor data noise suppression method and system based on mixed norm fitting - Google Patents
Aeromagnetic gradient tensor data noise suppression method and system based on mixed norm fitting Download PDFInfo
- Publication number
- CN109856690B CN109856690B CN201910150563.1A CN201910150563A CN109856690B CN 109856690 B CN109856690 B CN 109856690B CN 201910150563 A CN201910150563 A CN 201910150563A CN 109856690 B CN109856690 B CN 109856690B
- Authority
- CN
- China
- Prior art keywords
- compensation
- data
- change rate
- signal
- time change
- 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.)
- Expired - Fee Related
Links
Images
Abstract
The invention discloses a method and a system for suppressing noise of aeromagnetic gradient tensor data based on mixed norm fitting. The method comprises the following steps: establishing a magnetic gradient data comprehensive compensation model; obtaining a time change rate signal of observation data according to the magnetic gradient comprehensive compensation model; constructing a compensation equation based on the signal change rate according to the time change rate signal; establishing a data fitting objective function according to the compensation equation based on the signal change rate; solving the fitting objective function to obtain a compensation parameter; determining compensated magnetic gradient data according to the compensation parameters; and evaluating the compensation effect according to the compensated magnetic gradient data. The signal-to-noise ratio of the final compensation signal can be improved by adopting the method or the device.
Description
Technical Field
The invention relates to the technical field of geophysical aeromagnetic prospecting, in particular to a method and a system for suppressing noise of aeromagnetic gradient tensor data based on mixed norm fitting.
Background
Aeromagnetic exploration is an important geophysical exploration technique which is widely applied to mineral and environmental exploration. In recent years, magnetic gradient measurement systems based on superconducting quantum interferometers (SQUIDs) have been developed vigorously to enable direct measurement of high-precision magnetic gradient data, and superconducting full-tensor magnetic gradiometers are considered to be the development direction of third-generation aeromagnetic detection due to the advantages of high information content, high magnetic field sensitivity, small size and the like. At present, a full-tensor magnetic gradient measurement platform prototype based on SQUID is built in China, the full-tensor magnetic gradient measurement platform prototype is different from a traditional aeromagnetic gradient measurement flight platform which adopts an airborne mode to measure, a magnetic gradient measurement instrument and related equipment (GPS, an Inertial Navigation System (INS), a read-in read-out circuit and the like) are loaded in a nacelle by the prototype, and flight measurement is carried out in a helicopter dragging mode, the influence of a helicopter on magnetic measurement interference is greatly reduced by the prototype, even the influence of the helicopter can be ignored under the condition that a dragging cable is long enough, but a large amount of equipment which can generate magnetic interference still exists in the nacelle, so that the measurement precision of the flight platform is effectively improved through magnetic compensation treatment, and the high-resolution magnetic gradient measurement has important significance.
Because the magnetic gradient tensor measurement system adopts the high sensitivity characteristic of the SQUID sensor, the influence of the geomagnetic field gradient on the sensor in low-altitude exploration flight cannot be ignored, and the compensation coefficient obtained by the traditional high-altitude test flight cannot compensate the interference influence of the magnetic gradient signal in low-altitude flight on the sensor, so that the magnetic gradient processing data with the accuracy meeting the requirement cannot be obtained by directly using the traditional compensation method (namely, a mode of directly compensating the low-altitude exploration flight data by using the compensation coefficient obtained by using the high-altitude flight data). Therefore, in order to solve the problem, the original compensation operation mode needs to be changed, and the low-altitude flight survey data needs to be directly compensated. However, the important problem to be faced in this process is that the magnetic gradient signals inducing the interference source may also be effective signals for the survey purpose, and how to retain these effective signals on the premise of suppressing the interference and improve the signal-to-noise ratio of the final reconstructed signal becomes an important problem in the magnetic gradient data processing.
Disclosure of Invention
The invention aims to provide a method and a system for suppressing noise of aeromagnetic gradient tensor data based on mixed norm fitting, which can improve the signal-to-noise ratio of a final compensation signal.
In order to achieve the purpose, the invention provides the following scheme:
a noise suppression method for aeromagnetic gradient tensor data based on mixed norm fitting comprises the following steps:
establishing a magnetic gradient data comprehensive compensation model;
obtaining a time change rate signal of observation data according to the magnetic gradient comprehensive compensation model;
constructing a compensation equation based on the signal change rate according to the time change rate signal;
establishing a data fitting objective function according to the compensation equation based on the signal change rate;
solving the fitting objective function to obtain a compensation parameter;
determining compensated magnetic gradient data according to the compensation parameters;
and evaluating the compensation effect according to the compensated magnetic gradient data.
Optionally, obtaining a time rate of change signal of the observation data according to the magnetic gradient comprehensive compensation model specifically includes:
constructing a compensation equation Ax ═ b according to the magnetic gradient comprehensive compensation model;
obtaining time change rate signals DA and Db of the observed data according to the compensation equation;
wherein, a is directly acquired observation data, b is observation data for compensating a gradient value, x is a compensation parameter used in actual calculation, DA is a time change rate signal of the directly acquired observation data, and Db is a time change rate signal of the observation data for compensating a gradient value.
Optionally, the constructing a compensation equation based on the signal change rate according to the time change rate signal specifically includes:
according to the time change rate signal, constructing a compensation equation DAx ═ Db based on the signal change rate;
where x is a compensation parameter used in actual calculation, DA is a time change rate signal of observation data directly acquired, and Db is a time change rate signal of observation data for compensating a gradient value.
Optionally, the solving the fitted objective function to obtain the compensation parameter specifically includes:
acquiring initial compensation parameters and initial compensation data based on an Lp norm data fitting method;
calculating a first order time rate of change of the initial compensation data;
analyzing the regional characteristics of the effective signals according to the first-order time change rate to obtain an analysis result;
determining values of Lp norm orders p and q according to the analysis result;
and solving the data fitting objective function through a reweighted Gauss-Newton iterative algorithm according to the values of the p and the q to obtain compensation parameters.
Optionally, the evaluating the compensation effect according to the compensated magnetic gradient data specifically includes:
and according to the compensated magnetic gradient data, evaluating the compensation effect by adopting standardized parameters, wherein the standardized parameters comprise root mean square, standard deviation and improvement ratio.
A system for suppressing noise in aeromagnetic gradient tensor data based on mixed norm fitting, comprising:
the compensation model establishing module is used for establishing a magnetic gradient data comprehensive compensation model;
the time change rate signal acquisition module is used for acquiring a time change rate signal of observation data according to the magnetic gradient comprehensive compensation model;
the compensation equation determining module is used for constructing a compensation equation based on the signal change rate according to the time change rate signal;
the target function determining module is used for establishing a data fitting target function according to the compensation equation based on the signal change rate;
the solving module is used for solving the fitting objective function to obtain a compensation parameter;
the magnetic gradient data determining module is used for determining compensated magnetic gradient data according to the compensation parameters;
and the evaluation module is used for evaluating the compensation effect according to the compensated magnetic gradient data.
Optionally, the time change rate signal obtaining module specifically includes:
the first compensation equation determining unit is used for constructing a compensation equation Ax ═ b according to the magnetic gradient comprehensive compensation model;
the time change rate signal acquisition unit is used for obtaining time change rate signals DA and Db of the observation data according to the compensation equation;
wherein, a is directly acquired observation data, b is observation data for compensating a gradient value, x is a compensation parameter used in actual calculation, DA is a time change rate signal of the directly acquired observation data, and Db is a time change rate signal of the observation data for compensating a gradient value.
Optionally, the compensation equation determining module specifically includes:
the second compensation equation determining unit is used for constructing a compensation equation DAx ═ Db based on the signal change rate according to the time change rate signal;
where x is a compensation parameter used in actual calculation, DA is a time change rate signal of observation data directly acquired, and Db is a time change rate signal of observation data for compensating a gradient value.
Optionally, the solving module specifically includes:
acquiring initial compensation parameters and initial compensation data based on an Lp norm data fitting method;
calculating a first order time rate of change of the initial compensation data;
analyzing the regional characteristics of the effective signals according to the first-order time change rate to obtain an analysis result;
determining values of Lp norm orders p and q according to the analysis result;
and solving an objective function fitted by the self-adaptive Lp norm data through a re-weighted Gauss-Newton iterative algorithm according to the values of p and q.
Optionally, the evaluation module specifically includes:
and the evaluation unit is used for evaluating the compensation effect by adopting standardized parameters according to the compensated magnetic gradient data, wherein the standardized parameters comprise root mean square, standard deviation and improvement ratio.
According to the specific embodiment provided by the invention, the invention discloses the following technical effects: the invention provides a hybrid norm fitting-based aeromagnetic gradient tensor data noise suppression method, which can effectively improve the signal-to-noise ratio of a finally compensated magnetic gradient signal, so that a compensation algorithm can be conveniently and directly suitable for data compensation of a high-sensitivity magnetic gradient sensor.
Drawings
In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings needed to be used in the embodiments will be briefly described below, and it is obvious that the drawings in the following description are only some embodiments of the present invention, and it is obvious for those skilled in the art to obtain other drawings without inventive exercise.
FIG. 1 is a flow chart of a noise suppression method for aeromagnetic gradient tensor data based on mixed norm fitting according to the present invention;
FIG. 2 is a graph of magnetic susceptibility distribution values according to the present invention;
FIG. 3 is a comparison of the G1 channel measurement signals and different compensation data values of the magnetic gradient sensor of the present invention;
FIG. 4 is a comparison of the G2 channel measurement signals and different method compensation data values for the magnetic gradient sensor of the present invention;
FIG. 5 is a comparison of the G3 channel measurement signals and different compensation data values of the magnetic gradient sensor of the present invention;
FIG. 6 is a comparison of the G4 channel measurement signals and different method compensation data values for the magnetic gradient sensor of the present invention;
FIG. 7 is a comparison of the G5 channel measurement signals and different method compensation data values for the magnetic gradient sensor of the present invention;
FIG. 8 is a comparison of the G6 channel measurement signals and different method compensation data values for the magnetic gradient sensor of the present invention;
FIG. 9 is a block diagram of a noise suppression system for aeromagnetic gradient tensor data based on mixed norm fitting according to the present invention;
FIG. 10 is a graph illustrating the results of a compensation and filtering process on the magnetic gradiometer measurements of the actual G1 channel by the algorithm of the present invention, wherein (a) shows the comparison of pre-compensation and post-compensation data; (b) showing the comparison result of the filtered data before and after the compensation processing;
FIG. 11 is a Fourier frequency amplitude spectrum analysis comparison graph of the measured data, the compensated data and the filtered data after compensation and filtering processing are carried out on the magnetic gradiometer measured value of the G1 channel.
Detailed Description
The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.
The invention aims to provide a method and a system for suppressing noise of aeromagnetic gradient tensor data based on mixed norm fitting, which can improve the signal-to-noise ratio of a final compensation signal.
In order to make the aforementioned objects, features and advantages of the present invention comprehensible, embodiments accompanied with figures are described in further detail below.
Fig. 1 is a flowchart of a noise suppression method for aeromagnetic gradient tensor data based on mixed norm fitting. As shown in fig. 1, a noise suppression method for aeromagnetic gradient tensor data based on mixed norm fitting includes:
step 101: establishing a magnetic gradient data comprehensive compensation model;
according to a general assumption, the signals measured by the gradiometer consist of true signals, low frequency interference noise and high frequency electromagnetic noise, wherein the low frequency interference noise mainly originates from external magnetic source interference (residual magnetism, magnetic sensitivity, eddy current magnetic field of the measuring platform), and residual sensor system bias and parasitic magnetic field which are not calibrated, and the like. Therefore, when considering that the interference source of the measuring platform mainly comes from the interior of the nacelle and is rigidly connected with the magnetic sensor, the measurement gradient value of the magnetic sensor can be represented by the following model:
wherein l is a measurement channel number;
is the measured magnetic gradient value;
the true magnetic gradient value at the channel measurement point is obtained; heFor local background magnetic field, it can be approximated by global geomagnetic field model IGRF-12 model or local geomagnetic field modelTo; dHe＝(dHx,dHy,dHz)＝(He(t+1/fs)-He(t)). fs is a background magnetic field time rate of change signal, wherein fs is the data sampling frequency; b isrCompensating the corrected three-axis magnetometer measurements;
for the remaining systematic deviations that are not calibrated,
the residual magnetic interference of the external interference source is a constant value corresponding to the measured value of a single gradiometer within a certain time, and can be integrated into a deviation coefficient ol；
The linear operator is a linear operator of the magnetic induction interference of an external interference source and consists of (2p +1) multiplied by 3 coefficients;
the linear operator is the eddy current interference of an external interference source and consists of (2p +1) multiplied by 3 coefficients; h is0The linear operator is a parasitic magnetic field interference and consists of (2q +1) multiplied by 3 coefficients; n is high frequency Gauss white noise. By sorting the acquired signal model, a calculation model for magnetic gradient data compensation can be obtained, as follows:
wherein A is the amount of observation data and related calculation
Forming an NxM-level matrix, wherein N is the sampling point number of the gradiometer;
for the compensation parameters used in the actual calculation, a total of M is 12p +6q + 10. Because the model considers the system among various sensorsDelay characteristics (inertial navigation system, gradiometer and magnetometer), so that magnetic induction and eddy interference fields of external interference sources calculated by relating to different sensor measurement data and parasitic magnetic field interference are subjected to linear filtering approximation processing; when the system delay characteristic is not considered, p ═ q ═ 0 may be taken, and the model may be simplified to have M ═ 10 compensation coefficients x ═ o (o)l,ax,ay,az,bx,by,bz,cx,cz,cz)TThe calculation model of (2) can be expressed as:
step 102: obtaining a time change rate signal of observation data according to the magnetic gradient comprehensive compensation model;
constructing a compensation equation Ax ═ b according to the magnetic gradient comprehensive compensation model;
obtaining time change rate signals DA and Db of the observed data according to the compensation equation;
wherein, a is directly acquired observation data, b is observation data for compensating a gradient value, x is a compensation parameter used in actual calculation, DA is a time change rate signal of the directly acquired observation data, and Db is a time change rate signal of the observation data for compensating a gradient value.
Step 103: constructing a compensation equation based on the signal change rate according to the time change rate signal; the method specifically comprises the following steps:
according to the time change rate signal, constructing a compensation equation DAx ═ Db based on the signal change rate;
where x is a compensation parameter used in actual calculation, DA is a time change rate signal of observation data directly acquired, and Db is a time change rate signal of observation data for compensating a gradient value.
Step 104: establishing a data fitting objective function according to the compensation equation based on the signal change rate;
FIG. 2 is a graph of magnetic susceptibility distribution values according to the present invention. According to theoretical model calculation, if a target signal source is composed of a plurality of magnetic sources with different morphological characteristics, after a series of coordinate transformation (a geodetic coordinate system, a measuring platform coordinate system and a single gradiometer sensor measure parallel directions), magnetic anomaly and magnetic gradient anomaly signals have morphological similarity with the signals under the geodetic coordinate system, and only the middle part with larger amplitude of the magnetic anomaly has small amplitude fluctuation related to flight attitude change. The magnetic anomaly signal can thus be divided substantially into two regions, a first region being a relatively smooth flat region of the non-magnetic source region or the far magnetic source region, and a second region being a magnetic field fluctuation region of the magnetic source region or the magnetic source boundary region. In order to maintain a high signal-to-noise ratio of the compensated signal, it is desirable to minimize the ineffective fluctuations of the signal, especially the low frequency fluctuations in the signal, in the first region, since the high frequency fluctuations can be processed by low pass filtering without affecting the final compensation effect, while in the second region, it is desirable to maintain the amplitude and shape of the fluctuations of the signal, and the small signal oscillations associated with the change in flight attitude in the region, as much as possible.
In general, in a conventional L2 norm data fitting manner, because it is assumed that an error between a calculated value and an observed quantity satisfies Gaussian distribution, an interference signal can be well fitted when an effective signal amplitude is very low, and a final compensation signal tends to be uniform oscillation or a constant value with a small amplitude; however, when the calculated value and the residual quantity directly from the observed quantity have obvious amplitude variation (not satisfying Gaussian distribution), because the method can carry out indiscriminate fitting on the effective signal, the compensated data finally has obvious invalid signal oscillation near the effective signal area, and the signal-to-noise ratio of the compensation signal is reduced. And the Lp norm (p is more than or equal to 0 and less than or equal to 1) data fitting can effectively keep the step performance of the signal and reduce the invalid oscillation at the boundary of the reconstructed signal. In addition, the derivative information of the signal can further show the boundary characteristics of the signal, so that according to the analysis, the following objective function is constructed by combining the Lp norm and the time derivative information of the signal:
wherein p and q are the order of Lp norm,
or is
In actual calculation, the p and q values used by each region can be adaptively selected according to the characteristics of the signals, α is a fitting balance coefficient used for balancing the weight of two parts of fitting data, and can be generally selected
Since the objective function is not trivial in the case of p, q < 2, the discrete form of the objective function is usually approximated for ease of computation as a weighted L2 norm form of the objective function, namely:
where Φ (x) is the objective function, W1(x),W2(x) Is a norm weighting matrix whose diagonal elements are w1ii,w2iiThe specific definition mode is shown as the following formula:
and has anp,εqAll the gradient operator matrixes are constants with small enough values, and can take values from 1e to 8 in practical calculation, and the gradient operator matrixes used for calculating the time change are expressed as follows:
the above problem can be solved iteratively by a re-weighted gaussian-newtons-Newton method.
Step 105: solving the fitting objective function to obtain a compensation parameter;
step 105 specifically includes:
acquiring initial compensation parameters and initial compensation data based on an Lp norm data fitting method;
calculating a first order time rate of change of the initial compensation data;
analyzing the regional characteristics of the effective signals according to the first-order time change rate to obtain an analysis result;
determining values of Lp norm orders p and q according to the analysis result;
and solving the data fitting objective function through a reweighted Gauss-Newton iterative algorithm according to the values of the p and the q to obtain compensation parameters.
The specific steps of step 105 are:
step 1:
because the fitting objective function adopts different Lp fitting norm orders in different data areas for the measured data and the change rate of the measured data to improve the signal-to-noise ratio of the final compensation signal, an initial compensation data is required to be used as an algorithm input to judge the initial possibility area where the magnetic anomaly source is located. In the method, initial compensation data are given by using a result of general L2 norm fitting, a solution method can use a ridge regression algorithm or a solution method based on Principal Component Analysis (PCA), and a specific calculation scheme is as follows:
(1) the data fit objective function based on the L2 norm is:
(2) the corresponding compensation parameter x can be obtained by a ridge regression methodβAnd calculating to obtain compensation data
The basic calculation method comprises the following steps:
wherein β is a ridge regression parameter, and I is an M-order identity matrix.
(3) Through a solution method based on PCA, the calculation is more stable when the linear correlation of a system matrix is serious, and the specific calculation steps are as follows:
(i) the coefficient matrix A is subjected to mean value normalization processing column by column to obtain
Wherein, aiIs a column vector of the matrix a,
is a vector aiN is the number of sampling points.
(ii) To pair
Carrying out characteristic value decomposition on C ═ V ∑ VTObtaining the eigenvector matrix V ═ (V)1,v2,K,vM) And a corresponding eigenvalue λ1＞λ2＞K＞λM,(Σ＝diag(λ1,λ2,K,λM))
(iii) Selecting the eigenvalue to be greater than lambdaminM (M) components, forming a principal component transform compensation coefficient matrix P ═ AVm,Vm＝(v1,v2,K,vm) (ii) a Where P is the principal component matrix of A, Vm＝(v1,v2,K,vm) For a sub-matrix consisting of the first m eigenvector components, viFor the ith bit in the corresponding eigenvector matrix VAnd (5) sign vectors.
(iv) Solving principal component transformation compensation parameter y and corresponding compensated data G by least square methodrl0：
(II) step 2:
fitting the compensated data obtained from the correlation calculation to the initial L2 norm data
Calculating, judging and calibrating signal
The specific calculation method of the non-signal area, the smooth signal area and the non-smooth signal area in (1) is as follows:
(1) first, to the compensated signal
And rate of change thereof
Carrying out normalization processing to obtain normalized data di,bi：
(2) Setting a threshold theta of the normalized signal d, b1,θ2And judging and evaluating the order p, q of the fitting norm of the data and the change rate thereof according to the following rules:
this parameter will be used in the following data fitting calculation.
(III) step 3:
according to the order p and q of the fitting norm determined by the method, an objective function based on Lp mixed norm data fitting is solved through reweighting Gauss-Newton iteration, and a specific derivation process is shown as follows.
In order to find the minimum value of the objective function, it needs to satisfy the gradient value of the objective function equal to 0, due to the weighting term W in the original objective function1(x),W2(x) Related to the compensation parameter x, the approximate objective function of the original problem can be obtained through first-order Taylor expansion, and further, a corresponding iterative calculation formula can be obtained:
according to the condition that the minimum point of the objective function exists:
the following iterative calculation formula can be further derived:
the coefficient matrix and the right-end term in the linear iterative formula can be specifically expressed as:
in the formula
Is a norm weighting matrix whose diagonal elements are w1ii,w2iiThe specific definition mode is shown as the following formula:
and has anp,εqAll the constants are constants with small enough values, and the values can be 1e-8 in actual calculation.
Therefore, according to the above iterative calculation formula, the compensation parameter x and the corresponding compensated data can be calculated by the following algorithm:
(i) giving an initial value x0Iteration maximum step number maxstep and stagnation error;
(ii)doi＝1,...,maxstep；
(ii-1) updating the weighting coefficient matrix
(ii-2) constructing a Gauss-Hessian matrix HnAnd a gradient vector gn：
(ii-3) solving the equation set H by iterative linear equation set solving algorithm (e.g. PCG/GMRES/BiCG, etc.) or direct method (LU decomposition, etc.)nΔx＝-gnObtaining an updated variable Deltax and updating the calculated value xn+1/2＝xn+Δx；
(ii-4) performing a constrained interval projection x on the calculated updated variable according to the constraint conditionn+1＝Π|Q(xn+1/2)；
(ii-5) checking whether iteration stalls, if | | x is satisfiedn+1-xn||/||xnIf | | < error, iteration is stopped, otherwise, the step is shifted to (ii-1) to continue calculation.
(iii) Obtaining compensated data by calculating compensation parameter x
It should be noted that, in the calculation, the data amount and the data type of the coefficient matrix a used for constructing the magnetic compensation model equation are different (related to the modeling mode of the interference source model), the formed coefficient matrix a has serious ill-conditioned performance (i.e., the condition number is large), and in order to improve the stability of the solution parameters and the reusability of the compensation parameters, the principal component matrix P of the coefficient matrix a can be extracted by using the PCA analysis method in the above steps 1 and 3, and the data fitting objective function is constructed in the same way, and the compensation parameters y and the compensated data are obtained by solving
The specific objective function and compensation data solving steps are as follows:
step 106: determining compensated magnetic gradient data according to the compensation parameters;
step 107: and evaluating the compensation effect according to the compensated magnetic gradient data.
The compensation results are generally evaluated on the compensated data using the following normalization parameters: root mean square RMS, standard deviation σcImproving the specific IR. The specific definition of the parameter values is as follows:
wherein d iscTo compensate the processed data, d0N is the number of data of the same type for the original collected data.
In the case of a theoretical model test, since the exact signal value d is knowntTherefore, the residual r of the compensated signal and the accurate signal can be obtainedc＝dt-dcTo quantify the effectiveness of the algorithm, the following three evaluation parameters are specifically defined: root mean square RMS, standard deviation σ of residualcSignal-to-noise ratio SNR, and peak signal-to-noise ratio PSNR, the specific definitions of the parameter values are as follows:
wherein, the root mean square RMS of the residual quantity, the standard deviation sigma of the residual quantitycFor measuring parameter errors, the smaller the value is, the better the value is; the SNR, and PSNR, are used to measure the reconstruction accuracy of the effective signal, with larger values being better.
FIG. 9 is a block diagram of the noise suppression system for the aeromagnetic gradient tensor data based on the mixed norm fitting. As shown in fig. 9, a noise suppression system for the data of the aeromagnetic gradient tensor based on the mixed norm fitting includes:
a compensation model establishing module 201, configured to establish a magnetic gradient data comprehensive compensation model;
a time change rate signal obtaining module 202, configured to obtain a time change rate signal of the observation data according to the magnetic gradient comprehensive compensation model;
a compensation equation determining module 203, configured to construct a compensation equation based on the signal change rate according to the time change rate signal;
an objective function determining module 204, configured to establish a data fitting objective function according to the compensation equation based on the signal change rate;
a solving module 205, configured to solve the fit objective function to obtain a compensation parameter;
a magnetic gradient data determining module 206, configured to determine compensated magnetic gradient data according to the compensation parameter;
and the evaluation module 207 is used for evaluating the compensation effect according to the compensated magnetic gradient data.
The time change rate signal obtaining module 202 specifically includes:
the first compensation equation determining unit is used for constructing a compensation equation Ax ═ b according to the magnetic gradient comprehensive compensation model;
the time change rate signal acquisition unit is used for obtaining time change rate signals DA and Db of the observation data according to the compensation equation;
wherein, a is directly acquired observation data, b is observation data for compensating a gradient value, x is a compensation parameter used in actual calculation, DA is a time change rate signal of the directly acquired observation data, and Db is a time change rate signal of the observation data for compensating a gradient value.
The compensation equation determining module 203 specifically includes:
the second compensation equation determining unit is used for constructing a compensation equation DAx ═ Db based on the signal change rate according to the time change rate signal;
where x is a compensation parameter used in actual calculation, DA is a time change rate signal of observation data directly acquired, and Db is a time change rate signal of observation data for compensating a gradient value.
The solving module 205 specifically includes:
acquiring initial compensation parameters and initial compensation data based on an Lp norm data fitting method;
calculating a first order time rate of change of the initial compensation data;
analyzing the regional characteristics of the effective signals according to the first-order time change rate to obtain an analysis result;
and determining the values of Lp norm orders p and q according to the analysis result.
And solving an objective function fitted by the self-adaptive Lp norm data through a re-weighted Gauss-Newton iterative algorithm according to the values of p and q.
The evaluation module 207 specifically includes:
and the evaluation unit is used for evaluating the compensation effect by adopting standardized parameters according to the compensated magnetic gradient data, wherein the standardized parameters comprise root mean square, standard deviation and improvement ratio.
Specific example 1:
the effect of the magnetic gradient data compensation method of the present invention is verified and explained by theoretical model data and measured data. The SQUID magnetic gradient measuring instrument related in the invention has 6 magnetic gradient sensor output channels: G1-G6. In the actual measurement data verification, the present embodiment uses the data measured by the G1 channel to explain the magnetic gradient data compensation effect.
The effect of the magnetic gradient data compensation method proposed by the present invention is first verified and explained by means of simulation data. Assuming that the magnetic susceptibility value of the bedrock is 1e-8(SI), two magnetic source bodies with different magnetic susceptibility values and different forms are embedded in the bedrock, the left inclined plate-shaped body is an inhomogeneous source, the magnetic susceptibility values of the left inclined plate-shaped body are 1e-3, 2e-3 and 3e-3(SI) from top to bottom, and the magnetic susceptibility value of the right vertical tube-shaped body is 5e-3(SI), as shown in FIG. 3. Under a certain observation height, the magnetic anomaly signal coordinate is converted into a flight measurement coordinate system through attitude and position information actually measured by an inertial navigation system and a GPS, an approximate earth background magnetic field signal F0, a magnetic anomaly field three-component Fv and a magnetic gradient component signal Gi are obtained, a magnetic anomaly value is set to be Fv, and corresponding local magnetic field data can be expressed as Fa-F0 + Fv. Assuming that the interference source of the measuring platform is equivalent to a unit dipole, and by using the assumed remanence, magnetic induction and eddy current parameters and the unbalance coefficient, the actual measurement gradient value containing the aircraft maneuvering noise can be constructed by the following model construction:
wherein T is a coordinate transformation matrix from a geodetic coordinate system to a sensor measurement plane, and the used data for constructing the compensation model coefficient matrix comprises: f0, dF 0. The objective function solution method was fitted with three types of data: ridge regression, PCA analysis, the method of the present invention, obtained compensated magnetic gradient data are shown in FIGS. 5 and 6, where FIG. 3 shows the comparison of the G1 channel measurement signals of the magnetic gradient sensor of the present invention and the values of the compensated data of the different methods. FIG. 3 shows different curves respectively representing compensation data calculated by different methods, and FIG. 4 shows a comparison of the G2 channel measurement signal and compensation data values of different methods of the magnetic gradient sensor of the present invention; FIG. 5 is a comparison of the G3 channel measurement signals and different compensation data values of the magnetic gradient sensor of the present invention; FIG. 6 is a comparison of the G4 channel measurement signals and different method compensation data values for the magnetic gradient sensor of the present invention; FIG. 7 is a comparison of the G5 channel measurement signals and different method compensation data values for the magnetic gradient sensor of the present invention; FIG. 8 is a comparison of the G6 channel measurement signals and different method compensation data values for the magnetic gradient sensor of the present invention; the various evaluation indexes of the corresponding compensated data are shown in tables 1 to 3. A series of calculation results show that compared with the traditional method, the new method can effectively keep the amplitude and detail characteristics of the effective signal, and each evaluation parameter value is obviously superior to other two methods.
TABLE 1 data evaluation parameter statistics for magnetic gradient compensation method based on ridge regression data fitting solution
TABLE 2 data evaluation parameter statistics for magnetic gradient compensation method based on PCA analysis data fitting solution
TABLE 3 data evaluation parameter statistics for magnetic gradient compensation method based on new method data fitting solution
Then, the effect of the magnetic gradient data compensation method of the invention is verified and explained through the measured data. The sampled data height is 1km and therefore contains no effective magnetic gradient signals. FIG. 10 is a graph illustrating the results of a compensation and filtering process on the magnetic gradiometer measurements of the actual G1 channel by the algorithm of the present invention, wherein (a) shows the comparison of pre-compensation and post-compensation data; (b) and showing the comparison result between the filtered data and the filtered data after the compensation processing. FIG. 11 is a Fourier frequency amplitude spectrum analysis comparison graph of the measured data, the compensated data and the filtered data after compensation and filtering processing are carried out on the magnetic gradiometer measured value of the G1 channel. As can be seen from fig. 10 and 11, the magnetic compensation method can effectively suppress low-frequency interference in the measurement data, and the compensation data can further eliminate high-frequency noise therein after being filtered, so as to improve the compensation accuracy. Correspondingly, the quality evaluation indexes corresponding to the compensation and filtering processing results are given in table 4, and as can be seen by combining fig. 9-11 and table 4, the root mean square of the magnetic gradient data after the compensation and filtering processing is less than 20pT, and the improvement ratio IR can reach 2.0e3, so that a good compensation result can be obtained based on the method provided by the invention, and the requirement of instrument resolution can be met.
TABLE 4 data evaluation parameter statistics based on comprehensive magnetic gradient compensation method
So far, the basic method and the embodiment of the present invention have been described in detail with reference to the accompanying drawings. From the above description, those skilled in the art should clearly understand the magnetic gradient data compensation method based on single acquisition channel signals of SQUID aeromagnetic gradient meter according to the present invention. It is to be noted that, in the attached drawings or in the description, the implementation modes not shown or described are all the modes known by the ordinary skilled person in the field of technology, and are not described in detail.
The embodiments in the present description are described in a progressive manner, each embodiment focuses on differences from other embodiments, and the same and similar parts among the embodiments are referred to each other. For the system disclosed by the embodiment, the description is relatively simple because the system corresponds to the method disclosed by the embodiment, and the relevant points can be referred to the method part for description.
The principles and embodiments of the present invention have been described herein using specific examples, which are provided only to help understand the method and the core concept of the present invention; meanwhile, for a person skilled in the art, according to the idea of the present invention, the specific embodiments and the application range may be changed. In view of the above, the present disclosure should not be construed as limiting the invention.
Claims (6)
1. A noise suppression method for aeromagnetic gradient tensor data based on mixed norm fitting is characterized by comprising the following steps:
establishing a magnetic gradient data comprehensive compensation model;
obtaining a time change rate signal of observation data according to the magnetic gradient comprehensive compensation model;
constructing a compensation equation based on the signal change rate according to the time change rate signal;
establishing a data fitting objective function according to the compensation equation based on the signal change rate;
solving the fitted objective function to obtain compensation parameters, specifically comprising:
acquiring initial compensation parameters and initial compensation data based on an Lp norm data fitting method;
calculating a first order time rate of change of the initial compensation data;
analyzing the regional characteristics of the effective signals according to the first-order time change rate to obtain an analysis result;
determining values of fitting Lp norm orders p and q according to the analysis result;
solving the data fitting objective function through a reweighted Gauss-Newton iterative algorithm according to the values of p and q to obtain compensation parameters;
determining compensated magnetic gradient data according to the compensation parameters;
evaluating the compensation effect according to the compensated magnetic gradient data, which specifically comprises the following steps:
and according to the compensated magnetic gradient data, evaluating the compensation effect by adopting standardized parameters, wherein the standardized parameters comprise root mean square, standard deviation and improvement ratio.
2. The method for suppressing noise of the aeromagnetic gradient tensor data based on the mixed norm fitting as recited in claim 1, wherein the obtaining of the time change rate signal of the observed data according to the magnetic gradient comprehensive compensation model specifically comprises:
constructing a compensation equation Ax ═ b according to the magnetic gradient comprehensive compensation model;
obtaining time change rate signals DA and Db of the observed data according to the compensation equation;
wherein, a is directly acquired observation data, b is observation data for compensating a gradient value, x is a compensation parameter used in actual calculation, DA is a time change rate signal of the directly acquired observation data, and Db is a time change rate signal of the observation data for compensating a gradient value.
3. The method for suppressing noise of the aeromagnetic gradient tensor data based on the mixed norm fitting as recited in claim 1, wherein the constructing a compensation equation based on a signal change rate according to the time change rate signal specifically comprises:
according to the time change rate signal, constructing a compensation equation DAx ═ Db based on the signal change rate;
where x is a compensation parameter used in actual calculation, DA is a time change rate signal of observation data directly acquired, and Db is a time change rate signal of observation data for compensating a gradient value.
4. A noise suppression system for aeromagnetic gradient tensor data based on mixed norm fitting, comprising:
the compensation model establishing module is used for establishing a magnetic gradient data comprehensive compensation model;
the time change rate signal acquisition module is used for acquiring a time change rate signal of observation data according to the magnetic gradient comprehensive compensation model;
the compensation equation determining module is used for constructing a compensation equation based on the signal change rate according to the time change rate signal;
the target function determining module is used for establishing a data fitting target function according to the compensation equation based on the signal change rate;
the solving module is used for solving the fitting objective function to obtain a compensation parameter;
the magnetic gradient data determining module is used for determining compensated magnetic gradient data according to the compensation parameters;
the evaluation module is used for evaluating the compensation effect according to the compensated magnetic gradient data;
the solving module specifically includes:
acquiring initial compensation parameters and initial compensation data based on an Lp norm data fitting method;
calculating a first order time rate of change of the initial compensation data;
analyzing the regional characteristics of the effective signals according to the first-order time change rate to obtain an analysis result;
determining values of Lp norm orders p and q according to the analysis result;
according to the values of p and q, solving an objective function of self-adaptive Lp norm data fitting through a reweighted Gauss-Newton iterative algorithm;
the evaluation module specifically comprises:
and the evaluation unit is used for evaluating the compensation effect by adopting standardized parameters according to the compensated magnetic gradient data, wherein the standardized parameters comprise root mean square, standard deviation and improvement ratio.
5. The system for suppressing noise of magnetocaloric gradient tensor data based on mixed norm fitting as claimed in claim 4, wherein the time rate of change signal obtaining module specifically comprises:
the first compensation equation determining unit is used for constructing a compensation equation Ax ═ b according to the magnetic gradient comprehensive compensation model;
the time change rate signal acquisition unit is used for obtaining time change rate signals DA and Db of the observation data according to the compensation equation;
wherein, a is directly acquired observation data, b is observation data for compensating a gradient value, x is a compensation parameter used in actual calculation, DA is a time change rate signal of the directly acquired observation data, and Db is a time change rate signal of the observation data for compensating a gradient value.
6. The noise suppression system for the data of the aeromagnetic gradient tensor based on the mixed norm fitting as recited in claim 4, wherein the compensation equation determining module specifically comprises:
the second compensation equation determining unit is used for constructing a compensation equation DAx ═ Db based on the signal change rate according to the time change rate signal;
where x is a compensation parameter used in actual calculation, DA is a time change rate signal of observation data directly acquired, and Db is a time change rate signal of observation data for compensating a gradient value.
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 CN109856690A (en) | 2019-06-07 |
| CN109856690B true 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) |
Families Citing this family (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111079285B (en) * | 2019-12-16 | 2022-02-08 | 山东大学 | Full-tensor magnetic gradient data compensation optimization method and system |
| 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 |
| CN113534292B (en) * | 2021-07-09 | 2022-12-27 | 吉林大学 | Small signal model aeromagnetic compensation method based on forgetting factor RLS |
| CN115391726B (en) * | 2022-07-26 | 2023-10-27 | 睿励科学仪器（上海）有限公司 | Fitting optimization method and related device for ellipsometry system |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106959471B (en) * | 2017-04-21 | 2018-10-02 | 中国科学院电子学研究所 | Boat magnetic compensation method based on the non-linear boat total field gradient compensation model of magnetic |
| CN107291659B (en) * | 2017-05-16 | 2020-09-25 | 哈尔滨工程大学 | Recursive cosine transform method for extending plane modulus gradient field upwards in one step in plane geomagnetic abnormal field |
| 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
Also Published As
| Publication number | Publication date |
|---|---|
| CN109856690A (en) | 2019-06-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109856690B (en) | Aeromagnetic gradient tensor data noise suppression method and system based on mixed norm fitting | |
| CN109814163B (en) | Method and system for suppressing noise of aeromagnetic tensor data based on extended compensation model | |
| CN111433634B (en) | Magnetic compensation method based on aeromagnetic compensation error model | |
| CN109856689B (en) | Noise suppression processing method and system for superconducting aeromagnetic gradient tensor data | |
| CN107544042B (en) | Magnetometer array correction method | |
| Han et al. | A modified Tolles–Lawson model robust to the errors of the three-axis strapdown magnetometer | |
| Hardwick | Important design considerations for inboard airborne magnetic gradiometers | |
| Munschy et al. | Scalar, vector, tensor magnetic anomalies: Measurement or computation? | |
| CN106997035A (en) | A kind of gradometer bearing calibration based on magnetic gradient invariant | |
| CN110274586B (en) | Aeromagnetic compensation method including multi-optical system atom magnetometer direction error compensation | |
| CN110018429B (en) | Method and system for eliminating magnetic field interference caused by vibration of magnetic detection platform | |
| CN110187393B (en) | Aeromagnetic compensation method based on generalized regression neural network | |
| CN113156355B (en) | Magnetic interference compensation method of superconducting full tensor magnetic gradient measuring device | |
| CN113074721B (en) | Geomagnetic fingerprint construction method based on magnetic moment method | |
| CN108919156A (en) | Bearing calibration under three axis magnetometer line based on noise compensation | |
| CN109443337B (en) | Positioning navigation system and method based on NV color center in diamond | |
| Karshakov et al. | Aeromagnetic gradiometry and its application to navigation | |
| Ge et al. | Aeromagnetic system for a Multi-rotor unmanned aerial vehicle based on the Overhauser sensor | |
| Antti et al. | Separation of the geomagnetic variation field on the ground into external and internal parts using the spherical elementary current system method | |
| Wang et al. | Compensation for mobile carrier magnetic interference in a SQUID-based full-tensor magnetic gradiometer using the flower pollination algorithm | |
| Meyer et al. | SQUID technology for geophysical exploration | |
| Wan et al. | Improved component compensation for geomagnetic field vector measurement using Lagrange multiplier method | |
| CN114444299B (en) | Magnetic field reconstruction method based on distance weighted multipole expansion method | |
| CN108919368A (en) | A kind of system and method interfered for eliminating microsatellite remanent magnetism | |
| CN109655910A (en) | The two-parameter full waveform inversion method of Ground Penetrating Radar based on phasing |
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 |