CN107632964B - Downward continuation recursive cosine transform method for plane geomagnetic abnormal field - Google Patents
Downward continuation recursive cosine transform method for plane geomagnetic abnormal field Download PDFInfo
- Publication number
- CN107632964B CN107632964B CN201710795560.4A CN201710795560A CN107632964B CN 107632964 B CN107632964 B CN 107632964B CN 201710795560 A CN201710795560 A CN 201710795560A CN 107632964 B CN107632964 B CN 107632964B
- Authority
- CN
- China
- Prior art keywords
- formula
- plane
- abnormal field
- geomagnetic
- field
- 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.)
- Active
Links
Images
Landscapes
- Measuring Magnetic Variables (AREA)
- Complex Calculations (AREA)
Abstract
The invention discloses a downward continuation recursive cosine transform method for a plane geomagnetic abnormal field, and belongs to the field of geomagnetic field data processing. The method comprises the following specific steps: preprocessing the measurement data of the plane geomagnetic abnormal field; gridding the measurement data of the plane geomagnetic abnormal field; calculating a cosine transform spectrum of the planar geomagnetic abnormal field by using recursive cosine transform; setting a small amount epsilon according to the requirement of downward continuation precision, and selecting the value of a damping factor alpha; directly calculating the spectrum of the extended plane geomagnetic abnormal field by the cosine transform spectrum of the observation plane and the compensation factor of downward extension of the plane geomagnetic abnormal field; and performing one-dimensional recursive inverse cosine transformation on the wave spectrum in the row direction and the column direction to obtain a continuation value. The large-distance downward continuation of the geomagnetic abnormal field is stably carried out only by one-step spectrum transformation and inverse transformation transduction, so that the calculation amount of a downward continuation algorithm is reduced; the calculation amount of the recursive algorithm is small; the precision of large-distance downward continuation of the one-step wave number domain compensation of the plane geomagnetic abnormal field is improved.
Description
Technical Field
The invention belongs to the field of data processing of a geomagnetic field, and particularly relates to a downward continuation recursive cosine transform method for a planar geomagnetic abnormal field.
Background
In the field data processing, the downward continuation of geomagnetic anomaly has important significance for geological investigation, UXO detection, geomagnetic aided navigation and other aspects. In the areas with poor measurement conditions, such as the sea surface, the earth surface and the ultra-low altitude, the aviation geomagnetic anomaly measurement is indispensable. The aeromagnetic data is adopted to carry out downward continuation at different heights, so that the defect of magnetic field data in areas with severe conditions can be made up, the aeromagnetic data and ship survey geomagnetic anomaly data can be restored to the underwater by the downward continuation at a large distance, and a reference map is provided for the geomagnetic auxiliary navigation of an underwater vehicle. The downwardly extended bit field can highlight local geomagnetic anomaly, and improve the resolution of the anomaly and the reliability of data interpretation; downward continuation is also the core of many bit-field data processing techniques, such as the calculation of normalized overall gradients.
Downward continuation is a process of amplifying noisy data, and is an ill-defined problem, and if an inappropriate continuation method is adopted, a stable solution cannot be obtained. The solution obtained by the traditional downward continuation fast Fourier transform method is extremely unstable, and the reliable continuation depth generally does not exceed 2-3 times of the dot pitch. The downward continuation of the large-distance potential field becomes one of hot spots for processing the gravity magnetic data due to the important application value of the large-distance potential field. Dean studied the filtering characteristics of downward prolongation from the wavenumber domain, using a set of filter coefficients to match the frequency response of the downward prolongation operator (Dean W c. frequency analysis for gradient and magnetic interpretation, geophilics, 23(1): 97-127). Clarke applied wiener filter theory to design the downward continuation optimal filtering with interference (Clarke G KC. optimal second-derivative and downward-conversion filters, Geophysics,1969,34: 424-. Ku proposes correcting the downward-extending filter factor with a window function to eliminate the false anomaly caused by downward-extending amplified high-frequency interference (Ku C, Telford W M, Lim S H. the use of linear filtering in the diagrams, Geophysics,1971,36: 1174-. The characteristics of four regularized downward continuation filter factors were studied in Lima (regularization method of Lima. potential field downward continuation, geophysical reports, 1989,32 (5): 600-. Vertically projecting the actual measurement value of the bit field on the horizontal plane downwards to the continuation plane as an initial value of the continuation plane, and calculating the bit field value on the observation plane upwards by using fast Fourier transform; correcting the bit field value on the extension surface by using the difference between the measured value and the calculated value of the observation surface; the iteration is repeated until the difference value can be ignored (XushiZhe. comparison of the iterative method and the FFT method bit field downward continuation effect, the geophysical report, 2007, 50 (1): 285-. The Liudong Jia and the like research the spatial iterative method of Xueshi Zhe in the wave number domain, form a wave number domain iterative method and analyze the filtering characteristics, and strictly prove the convergence of the algorithm (Liudong Jia, Shangtian, Jiashihai and the like, the wave number domain iterative method of extending the bit field downwards and the convergence thereof, the geophysical report 2009, 52 (6): 1599-1605). Calf et al introduced the concept of cut-off wave numbers, less than which do downward continuation of the normal continuation factor, greater than which do downward continuation of the regularized continuation operator (Xiao Niu. an improved regularized downward continuation of positional field data, Journal of Applied geomatics, 2014,106: 114-. Mannational celebration et al use an iterative process of upward continuation and horizontal derivative combination to achieve stable downward continuation of bit field data (Guoqing Ma, Cai Liu, Danian Huang, Lili Li. A stable iterative downward continuation of positional field data, Journal of Applied geomatics, 2013,98: 205-reservoir 211). The method combines a wavenumber domain generalized inverse algorithm and an iterative method, and realizes the stable downward continuation of a potential field by a successive compensation method, wherein the continuation has higher precision (research and application of a compensation downward continuation method of the Gaoyu, Luoyao, Wenwu. the research and application of the geophysical science report 2012, 55(8):2747 + 2756).
Although the iterative calculation method and the compensation iteration method can perform stable large-distance downward continuation to obtain higher continuation accuracy, repeated iteration is needed, the calculation efficiency of the downward continuation algorithm is seriously influenced by the increase of the iteration times, and the ideal downward continuation can be approached under the condition of more iteration times. If the iteration times are too few, the accuracy and the robustness of the downward continuation algorithm are both seriously reduced.
The invention provides a large-distance downward extension recursive cosine transform method for compensating a one-step wave number domain of a plane geomagnetic abnormal field. And deducing a one-step wave number domain compensation operator of downward continuation based on cosine transform according to a wave number domain iterative compensation method of downward continuation. And (3) directly calculating the geomagnetic abnormal field of the continuation plane by using the recursive cosine transform and the inverse transform and by using the downward continuation one-step wave number domain compensation operator and the geomagnetic abnormal field of the observation plane without an iteration process. Compared with a downward continuation Fourier transform method of wave number domain iterative compensation, the downward continuation cosine transform method of one-step wave number domain compensation does not need repeated iteration; cosine transform has higher energy compression performance than fourier transform, and Gibbs boundary effect can be reduced. Therefore, compared with a downward continuation Fourier transform method of wave number domain iterative compensation and the like, the downward continuation recursive cosine transform method of one-step wave number domain compensation of the geomagnetic abnormal field not only can realize large-distance downward continuation, but also reduces the operation time and the calculation error of a continuation algorithm, and simplifies the calculation flow. The recursive cosine transform and the inverse transform can also adapt to data sequences with any length, and the fast fourier transform and the inverse transform generally require that the length of the data sequence is an integer power of a base number, and zero padding can increase the calculation amount of the algorithm.
Disclosure of Invention
The invention aims to provide a downward continuation recursive cosine transform method of a planar geomagnetic abnormal field, which improves the algorithm precision and can be applied to data sequences with any length.
The purpose of the invention is realized by the following technical scheme:
a downward continuation recursive cosine transform method for a plane geomagnetic abnormal field comprises the steps of preprocessing measurement data of the observation plane geomagnetic abnormal field, and removing measurement data with large errors; and meshing the geomagnetic abnormal field measurement data which are irregularly and discretely distributed on the observation plane by using a radial basis function interpolation method to obtain the meshed geomagnetic abnormal field measurement data. Then, converting the meshed geomagnetic anomaly measurement data into a cosine transform spectrum by using recursive cosine transform; utilizing a one-step wave number domain compensation operator for downward continuation to obtain a geomagnetic abnormal spectrum on a continuation plane; and then, calculating by using the inverse recursive cosine transform to obtain the geomagnetic abnormal field on the continuation plane. The specific steps are described below, where steps 5, 6 and 7 are possible for parallel computing.
Step 1, selecting gridding parameters, preprocessing the measurement data of the planar geomagnetic abnormal field, and removing coarse measurement errors which exceed a reasonable range.
Step 2, carrying out network on discrete geomagnetic abnormal field measurement data on an observation plane by using a radial basis function interpolation methodGridding to obtain gridded geomagnetic anomaly measurement data delta T (n) on observation planex,ny,z0)，nx＝1,2,…,Mx，ny＝1,2，…,My. The equation of the observation plane is z ═ z0The equation for the downward continuation plane is z ═ z0The + Δ z, z axis is vertically downward.
Step 3, one-dimensional recursive cosine transform of line direction and column direction is carried out on the geomagnetic abnormal field on the observation plane in sequence, the calculation orders of the two directions can be exchanged, and the two-dimensional cosine transform wave spectrum delta T of the geomagnetic abnormal field is obtainedC(ux,uy,z0)。
Step 4, setting a small quantity epsilon for representing the accuracy of the downward continuation algorithm, selecting the value of a damping factor alpha, wherein alpha is more than 0 and less than 1, usually taking alpha as 0.1, and solving the iteration number n in the iterative compensation method according to an inequality shown in formula (1)I。
In the formula
Step 5, according to the formula (2), based on the geomagnetic abnormal field spectrum Delta T of the observation planeC(ux,uy,z0) Calculating the spectrum Delta T of the geomagnetic abnormal field of the continuation planeC(ux,uy,z0+Δz)。
In the formula (I), the compound is shown in the specification,
a step wave number domain compensation factor for downward continuation of the plane geomagnetic abnormal field, which is expressed as
One-step wave number domain compensation factor is equivalent to n in wave number domain iterative compensation methodIAnd (5) performing secondary iteration compensation calculation.
Step 6, successively comparing delta TC(ux,uy,z0+ Δ z) is subjected to one-dimensional recursive inverse cosine transform in the row direction and the column direction, the calculation order of the two directions can be exchanged, and the continuation value Δ T (x, y, z) of the geomagnetic abnormal field on the continuation plane is obtained0+Δz)。
The invention has the beneficial effects that:
the large-distance downward extension recursive cosine transform method for compensating the one-step wave number domain of the plane geomagnetic abnormal field has the advantages of no need of iteration, small calculated amount, large extension distance, high extension precision and the like, and solves the problem that the compensation iteration large-distance downward extension method needs repeated iteration; the method only needs one step of spectrum transformation and inverse transformation transduction to stably carry out large-distance downward extension on the geomagnetic abnormal field, simplifies the algorithm steps of large-distance downward extension, and reduces the calculation amount of a downward extension algorithm; the bit field spectrum transformation and the inverse transformation adopt a recursive algorithm of cosine transformation and inverse cosine transformation, and the calculated amount is small; compared with fast Fourier transform, the Gibbs boundary effect is reduced, and the large-distance downward continuation precision of the one-step wave number domain compensation of the planar geomagnetic abnormal field is improved. The recursive cosine transform and the inverse transform can adapt to a data sequence with any length, and the fast fourier transform and the inverse transform generally require that the length of the data sequence is an integer power of a base number, and zero padding increases the calculation amount of the algorithm.
Drawings
FIG. 1 is a block diagram of a planar geomagnetic anomaly field downward continuation recursive cosine transform method
FIG. 2 is a flowchart of a downward continuation recursive cosine transform method for a planar geomagnetic anomaly field
FIG. 3 is a sphere model experiment result of large distance downward continuation based on one-step wave number domain compensation of Fourier transform and cosine transform
FIG. 4 is the experimental result of the mixed model of sphere and cuboid with large distance downward continuation based on one-step wave number domain compensation of Fourier transform and cosine transform
FIG. 5 is a large-distance downward continuation mean square error curve of ball and cuboid mixed model data compensated in one-step wave number domain
FIG. 6 is a relative error curve of large-distance downward continuation of ball-cuboid hybrid model data compensated in one-step wave number domain
Detailed Description
The following further describes embodiments of the present invention with reference to the accompanying drawings:
step 1, preprocessing the measurement data of the plane geomagnetic abnormal field. And setting a reasonable range of the measurement error and the geomagnetic abnormal value, and removing the measurement data of the gross error by regarding the geomagnetic abnormal field measurement data beyond the reasonable range as the gross error.
And 2, gridding the measurement data of the plane geomagnetic abnormal field. And selecting gridding parameters, and gridding the discrete geomagnetic abnormal field measurement data by using a radial basis function interpolation method to obtain the gridded geomagnetic abnormal field measurement data on the observation plane.
The radial basis function interpolation method is an interpolation method with high accuracy, and the basis function is utilized to determine the optimal weight from a known data point to an interpolation grid node; the method is suitable for interpolation calculation of a large amount of point data, has high prediction precision, and can better reflect the change of the data. The planar interpolation expression of the radial basis function of the point to be interpolated is
Wherein (x, y) is the coordinate of the interpolation point, f0(x,y)＝c1xx+c1yy+c0，(xj,yj) Is the coordinate of the sampling point, λjFor the weight corresponding to each sample point,
for radial basis functions, | | · | | is the euclidean norm. Selecting thin plate spline function asRadial basis function, d·j＝||x-xj,y-yjIf there is
The sampling point is given the equation constraint condition of
Assuming that f (x, y) has a second continuous derivative, the energy function is expressed as:
minimizing the energy function to obtain the required orthogonal condition of
To solve the coefficient c in the formula (1)0、c1x、c1yAnd the weight value lambdajThe combination of formula (3) and formula (5) yields a coefficient of c0、c1x、c1yAnd the weight value lambdajHas a linear equation set of
Solving a linear equation set with the coefficient matrix being a symmetric positive definite matrix to obtain a coefficient c0、c1x、c1yAnd the weight value lambdaj(j is 1,2, …, p), and the gridded geomagnetic abnormal field data on the observation surface is obtained by substituting the formula (1), nx＝1,2,…,Mx，ny＝1,2,…,MyThe equation of the observation plane is z ═ z0The z-axis is vertically downward.
Step 3, utilizing recursion remainderChord transformation by Δ T (n)x,ny,z0) Calculating cosine transform spectrum delta T of plane geomagnetic abnormal fieldC(ux,uy,z0). Definition of
Step 1) pressing formula (8) with Δ T (u)x,uy,z0) Computing
In the formula ux＝nx，θux＝(ux-1)π/Mx，
And
determined by the recursive formula (9).
In the formula, mx＝1,2,…,Kx，Kx＝[(Mx+1)/2]Symbol [. C]Indicating rounding.
Step 2) pressing formula (11) consisting of
Calculating Δ TC(ux,uy,z0)，
In the formula uy＝ny，θuy＝(uy-1)π/My，
And
determined by the recursive formula (12).
In the formula, my＝1,2,…,Ky，Ky＝[(My+1)/2]Symbol [. C]Indicating rounding.
And 4, setting a small amount epsilon according to the requirement of downward continuation precision, and selecting the value of a damping factor alpha, wherein alpha is more than 0 and less than 1, and is usually 0.1. When the iteration number is too small or too large, the downward continuation precision does not meet the requirement generally, and n satisfying the inequality shown in the formula (14) is solvedI。
In the formula
Step 1) k is equal to 0, and n is determined according to practical experienceIInitial interval of (2)
Step 2) selecting a section
Intermediate point of (2)
Computing
If it satisfies
Stopping searching; if not, then,
step 3) selecting the interval
Intermediate point of (2)
Computing
If it satisfies
Stopping searching; if not, then,
and 4) turning to step 2) to continue searching when k is equal to k + 1.
Step 5, according to the formula (15), from the cosine transform spectrum Delta T of the geomagnetic abnormal field on the observation surfaceC(ux,uy,z0) Directly calculating the wave spectrum delta T of the extended plane geomagnetic abnormal field by using the compensation factor of the one-step wave number domain extended downwards by the plane geomagnetic abnormal fieldC(ux,uy,z0+Δz)。
In the formula (I), the compound is shown in the specification,
a step wave number domain compensation factor for downward continuation of the plane geomagnetic abnormal field, which is expressed as
Step 6, for delta TC(ux,uy,z0- Δ z) of the geomagnetic anomaly field on the extended plane to obtain a linear vector of the geomagnetic anomaly field on the extended planex,ny,z0+Δz)。
Step 1) according to formula (17) by Δ TC(ux,uy,z0+ Δ z) calculation
In the formula (I), the compound is shown in the specification,
and
determined by recursive formula (18).
Step 2) pressing formula (19) consisting of
Calculating Δ T (n)x,ny,z0+Δz)，
In the formula (I), the compound is shown in the specification,
and
determined by the recursive formula (20).
And (3) utilizing recursive cosine transformation and inverse transformation to continue downwards from one-step wave number domain of the geomagnetic abnormal field of the observation plane by a large distance to obtain the geomagnetic abnormal field of the continuation plane. A block diagram of a cosine transform method for extending a planar geomagnetic abnormal field downwards by a large distance in a one-step wave number domain is shown in fig. 1. Fig. 2 shows a flowchart of a cosine transform method for extending downward a planar geomagnetic abnormal field in a one-step wave number domain by a large distance.
Two error indexes of a root mean square error and a relative error are defined for representing the extension error of the plane geomagnetic abnormal field extending downwards at a large distance in one-step wave number domain.
Root Mean Square Error (RMSE) for downward continuation using Fourier transformFIs composed of
Root Mean Square Error (RMSE) for downward continuation using cosine transformCIs composed of
Relative error RE of downward continuation using Fourier transformFIs composed of
Relative error RE of downward continuation using cosine transformCIs composed of
In the formula,. DELTA.T (n)x,ny,z0+ Δ z) is the model theoretical value of the extended planar geomagnetic anomaly field, Δ TF(nx,ny,z0+ Δ z) is an extension value of the extended planar geomagnetic abnormal field obtained using fourier transform and inverse transform.
Experiment one: z is the observation plane 00, the x and y coordinate ranges of the observation region are both [ -10000m,10000m]The number of meshes is 256 × 256, and the depth of downward continuation is Δ z 1568.6m (20 times dot pitch). The positions of the six spherical magnetic sources are respectively [0m,0m and 4000m]、[2000m,800m,3400m]、[1800m,-1600m,3500m]、[-2000m,-1800m,3600m]、[-1000m,2000m,3700m]And [1500m,0m,4300m]The magnetization intensity is 1.6 x 10 respectively5A/m、1.4×105A/m、1.3×105A/m、1.5×105A/m、1.2×105A/m and 1.8X 105A/m, the spherical radii are respectively 10m, 12m, 8m, 15m, 20m and 18m, the magnetic declination angles are respectively 45 degrees, 60 degrees, 75 degrees, 45 degrees, 60 degrees and 25 degrees, and the magnetic declination angles are respectively 45 degrees, 35 degrees, 15 degrees, 35 degrees, 30 degrees, 45 degrees and 45 degrees. The magnetic dip angle and the magnetic declination angle of the geomagnetic field are 60 degrees and 10 degrees respectively. The magnetic measurement noise on the observation surface is high-frequency-band Gaussian noise with zero mean and 20nT standard deviation.
A variation curve of the extended root mean square error of the downward extension of the sphere model data with the number of iterations is shown in fig. 3, where a black line is a root mean square error curve obtained by fourier transform and a red line is a root mean square error curve obtained by cosine transform. A curve of a downward continuation relative error with the number of iterations is shown in fig. 4, where a black line is a relative error curve obtained by fourier transform and a red line is a relative error curve obtained by cosine transform. The result shows that the continuation precision based on cosine transform and inverse transform is higher than that based on Fourier transform and inverse transform, and the precision of one-step iterative compensation is higher than that of non-iterative compensation; the error of the one-step iteration compensation algorithm is reduced along with the increase of the iteration times, is not reduced along with the increase of the iteration times, and is obviously increased for the relative error.
Experiment two: and adding two cuboid magnetic sources in the simulation experiment area. The first cuboid has central position of [0m,0m,3000m ], length, width and height of 80m, 40m and 40m, total magnetization of 1A/m, and declination angle of 50 deg. and 30 deg.. The second rectangular solid has a central position of [1000m,0m,3500m ], a length, a width and a height of 100m, 60m and 20m, respectively, a total magnetization of 0.8A/m, and a declination angle of 60 DEG and 15 deg, respectively.
A variation curve of the extended root mean square error of the model data of the mixture of the sphere and the rectangular parallelepiped, which is extended downward, with the number of iterations is shown in fig. 5, in which a black line is a root mean square error curve obtained by fourier transform, and a red line is a root mean square error curve obtained by cosine transform. A curve of a downward continuation relative error with the number of iterations is shown in fig. 6, where a black line is a relative error curve obtained by fourier transform and a red line is a relative error curve obtained by cosine transform. The result also shows that the continuation precision based on cosine transform and inverse transform is higher than that based on Fourier transform and inverse transform, and the precision of one-step iterative compensation is higher than that of non-iterative compensation; the error of the one-step iteration compensation algorithm is reduced along with the increase of the iteration number at first, and then is not obviously reduced along with the increase of the iteration number, and the relative error is also obviously increased.
The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention, and various modifications and changes may be made by those skilled in the art. Any modification, equivalent replacement, or improvement made within the spirit and principle of the present invention should be included in the protection scope of the present invention.
Claims (1)
1. A downward continuation recursive cosine transform method for a plane geomagnetic abnormal field is characterized by comprising the following steps
(1) Preprocessing the measurement data of the plane geomagnetic abnormal field, selecting gridding parameters, preprocessing the measurement data of the plane geomagnetic abnormal field, and removing coarse measurement errors which exceed a reasonable range;
(2) gridding the measurement data of the plane geomagnetic abnormal field by using a radial basis function interpolation method to gridde the measurement data of the discrete geomagnetic abnormal field on the observation plane to obtain the gridded geomagnetic abnormal measurement data delta T (n) on the observation planex,ny,z0)，nx＝1,2,…,Mx，ny＝1,2，…,My；
(3) Using a recursive cosine transform consisting of Δ T (n)x,ny,z0) Calculating cosine transform spectrum delta T of plane geomagnetic abnormal fieldC(ux,uy,z0)；
(4) Setting a small amount epsilon according to the requirement of downward continuation precision, selecting the value of a damping factor alpha, wherein alpha is more than 0 and less than 1, and taking alpha as 0.1;
(5) cosine transform spectrum Delta T of geomagnetic abnormal field from observation surfaceC(ux,uy,z0) Directly calculating the wave spectrum delta T of the extended plane geomagnetic abnormal field by using the compensation factor of the one-step wave number domain extended downwards by the plane geomagnetic abnormal fieldC(ux,uy,z0+Δz)；
(6) For Δ TC(ux,uy,z0+ Δ z) to obtain an extension value Δ T (n) of the geomagnetic abnormal field on the extension planex,ny,z0+Δz)；
The step (2) is specifically as follows:
the planar interpolation expression of the radial basis function of the point to be interpolated is
Wherein (x, y) is the coordinate of the interpolation point, f0(x,y)＝c1xx+c1yy+c0，(xj,yj) Is the coordinate of the sampling point, λjFor the weight corresponding to each sample point,
is a radial basis function, and | is | · | | | is a euclidean norm;
selecting thin plate spline function as radial basis function, let d·j＝||x-xj,y-yjIf there is
The sampling point is given the equation constraint condition of
Assuming that f (x, y) has a second continuous derivative, the energy function is expressed as:
minimizing the energy function to obtain the required orthogonal condition of
To solve the coefficient c in the formula (1)0、c1x、c1yAnd the weight value lambdajThe combination of formula (3) and formula (5) yields a coefficient of c0、c1x、c1yAnd the weight value lambdajHas a linear equation set of
Solving a linear equation set with the coefficient matrix being a symmetric positive definite matrix to obtain a coefficient c0、c1x、c1yAnd the weight value lambdaj(j is 1,2, …, p), and the gridded geomagnetic abnormal field data on the observation surface is obtained by substituting the formula (1), nx＝1,2,…,Mx，ny＝1,2,…,MyThe equation of the observation plane is z ═ z0The z-axis is vertically downward;
the step (3) is specifically as follows:
definition of
(3.1) by the formula (8) of Δ T (u)x,uy,z0) Calculating Δ T1 C(ux,ny,z0)，
In the formula ux＝nx，θux＝(ux-1)π/Mx，
And
determined by the recursive formula (9);
in the formula, mx＝1,2,…,Kx，Kx＝[(Mx+1)/2]Symbol [. C]Representing rounding;
(3.2) by the formula (11) of1 C(ux,ny,z0) Calculating Δ TC(ux,uy,z0)，
In the formula uy＝ny，θuy＝(uy-1)π/My，
And
determined by the recursive formula (12);
in the formula, my＝1,2,…,Ky，Ky＝[(My+1)/2]Symbol [. C]Representing rounding;
the step (4) is specifically as follows:
when the iteration number is too small or too large, the downward continuation precision does not meet the requirement, and n meeting the inequality shown in the formula (14) is solvedI
In the formula
(4.1) k is 0, and n is determined empiricallyIInitial interval of (2)
(4.2) selection of intervals
Intermediate point of (2)
Computing
If it satisfies
Stopping searching; if not, then,
(4.3) selection intervals
Intermediate point of (2)
Computing
If it satisfies
Then stop searching(ii) a If not, then,
(4.4) k is k +1, and turning to the step 2) to continue searching;
the step (5) is specifically as follows:
in the formula (I), the compound is shown in the specification,
a step wave number domain compensation factor for downward continuation of the plane geomagnetic abnormal field, which is expressed as
The step (6) is specifically as follows:
(6.1) by the formula (17) ofC(ux,uy,z0+ Δ z) calculating Δ T1 C(nx,uy,z0+Δz)，
In the formula (I), the compound is shown in the specification,
and
determined by recursive formula (18);
(6.2) by the formula (19) of1 C(nx,uy,z0+ Δ z) calculation Δ T (n)x,ny,z0+Δz)，
In the formula (I), the compound is shown in the specification,
and
determined by a recursive formula (20);
root Mean Square Error (RMSE) for downward continuation using Fourier transformFIs composed of
Root Mean Square Error (RMSE) for downward continuation using cosine transformCIs composed of
Relative error RE of downward continuation using Fourier transformFIs composed of
Relative error RE of downward continuation using cosine transformCIs composed of
In the formula,. DELTA.T (n)x,ny,z0+ Δ z) is the model theoretical value of the extended planar geomagnetic anomaly field, Δ TF(nx,ny,z0+ Δ z) is an extension value of the extended planar geomagnetic abnormal field obtained using fourier transform and inverse transform.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710795560.4A CN107632964B (en) | 2017-09-06 | 2017-09-06 | Downward continuation recursive cosine transform method for plane geomagnetic abnormal field |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710795560.4A CN107632964B (en) | 2017-09-06 | 2017-09-06 | Downward continuation recursive cosine transform method for plane geomagnetic abnormal field |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107632964A CN107632964A (en) | 2018-01-26 |
| CN107632964B true CN107632964B (en) | 2021-06-11 |
Family
ID=61100103
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710795560.4A Active CN107632964B (en) | 2017-09-06 | 2017-09-06 | Downward continuation recursive cosine transform method for plane geomagnetic abnormal field |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107632964B (en) |
Families Citing this family (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109190082B (en) * | 2018-08-15 | 2022-12-30 | 中国人民解放军61540部队 | Method for downwardly extending aviation gravity data |
| CN110389391B (en) * | 2019-08-01 | 2020-12-15 | 自然资源部第二海洋研究所 | Heavy magnetic bit field analytic extension method based on spatial domain |
| CN110543611B (en) * | 2019-08-15 | 2022-11-25 | 桂林理工大学 | Low latitude magnetic abnormal data magnetization pole calculation method and device |
| CN110531430A (en) * | 2019-08-29 | 2019-12-03 | 中国石油天然气集团公司 | Processing method, device and the electronic equipment of submarine pipeline magnetic survey data |
| CN111337992B (en) * | 2020-03-23 | 2021-04-06 | 兰州大学 | Method for obtaining depth of field source based on downward continuation of bit field data |
| CN114048636B (en) * | 2022-01-11 | 2022-05-03 | 中国测绘科学研究院 | Gravity anomaly calculation method and device based on wavelet transformation |
| CN116774301B (en) * | 2023-06-26 | 2024-05-14 | 中国人民解放***箭军工程大学 | Method, system, electronic equipment and medium for downward continuation of heavy magnetic bit field regularization |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103900564A (en) * | 2014-03-03 | 2014-07-02 | 哈尔滨工程大学 | Submergence assisted geomagnetic anomaly inversion velocity measurement/underwater continuous positioning method |
| CN104215259A (en) * | 2014-08-22 | 2014-12-17 | 哈尔滨工程大学 | Inertial navigation error correction method based on geomagnetism modulus gradient and particle filter |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP4656024B2 (en) * | 2006-08-22 | 2011-03-23 | 株式会社デンソー | Abnormality detection device for rotation angle detection device |
-
2017
- 2017-09-06 CN CN201710795560.4A patent/CN107632964B/en active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103900564A (en) * | 2014-03-03 | 2014-07-02 | 哈尔滨工程大学 | Submergence assisted geomagnetic anomaly inversion velocity measurement/underwater continuous positioning method |
| CN104215259A (en) * | 2014-08-22 | 2014-12-17 | 哈尔滨工程大学 | Inertial navigation error correction method based on geomagnetism modulus gradient and particle filter |
Non-Patent Citations (6)
| Title |
|---|
| 任意长离散余弦变换的快速递归算法;沈宏君;《电子与信息学报》;20070228;第29卷(第2期);第418-420页 * |
| 位场向下延拓三种迭代方法之比较;曾小牛 等;《地球物理学进展》;20110630;第26卷(第3期);第908-915页 * |
| 利用余弦变换计算重力异常的向上延拓;张凤旭 等;《地球物理学进展》;20070228;第22卷(第1期);第57-62页 * |
| 地磁异常数据向下延拓算法研究;马涛;《中国优秀硕士学位论文全文数据库信息科技辑》;20120715(第7期);第I138-20页 * |
| 基于离散余弦变换的位场谱方法及应用;张凤琴;《中国博士学位论文全文数据库基础科学辑》;20071015(第4期);正文第29-35页 * |
| 正则化薄板样条函数拟合地层界面;王宝龙 等，;《第十五届全国数学地质与地学信息学术研讨会》;20161021;第160-165页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107632964A (en) | 2018-01-26 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107632964B (en) | Downward continuation recursive cosine transform method for plane geomagnetic abnormal field | |
| CN107291659B (en) | Recursive cosine transform method for extending plane modulus gradient field upwards in one step in plane geomagnetic abnormal field | |
| CN110389391B (en) | Heavy magnetic bit field analytic extension method based on spatial domain | |
| CN110133644B (en) | Ground penetrating radar three-dimensional forward modeling method based on interpolation scale function method | |
| CN105319589B (en) | A kind of fully automatic stereo chromatography conversion method using local lineups slope | |
| CN105652320B (en) | Reverse-time migration imaging method and device | |
| CN113255230B (en) | Gravity model forward modeling method and system based on MQ radial basis function | |
| CN111856596A (en) | Layered medium resistivity anisotropy ocean controllable source electromagnetic rapid inversion method | |
| CN115292973A (en) | Arbitrarily sampled space wave number domain three-dimensional magnetic field numerical simulation method and system | |
| CN111999778A (en) | Antarctic continental Mohuo surface depth inversion method based on satellite gravity gradient data | |
| CN110413939B (en) | Arrival angle estimation method based on atomic norm | |
| CN107942374A (en) | Diffracted wave field extracting method and device | |
| CN113419280B (en) | Pre-stack crack density estimation method based on improved ellipse fitting | |
| CN114167511A (en) | Continuous-fraction expansion downward continuation-based bit field data rapid inversion method | |
| CN115508898B (en) | G-S conversion grounding long-wire source transient electromagnetic rapid forward and backward method and system | |
| CN106842297A (en) | Borehole restraint unstable state method for correcting phase | |
| CN109557588B (en) | Coal mine underground two-dimensional mine seismic wave velocity inversion dimension reduction method | |
| CN108802819B (en) | A kind of trapezoidal grid finite difference Simulation of Seismic Wave method of uniform depth sampling | |
| US20170108618A1 (en) | Geophysical inversion using sparse modeling | |
| Zhou et al. | Optimizing orthogonal-octahedron finite-difference scheme for 3D acoustic wave modeling by combination of Taylor-series expansion and Remez exchange method | |
| CN109188522B (en) | Velocity field construction method and device | |
| CN109387872B (en) | Surface multiple prediction method | |
| Morgan et al. | Inversion of combined gravity and bathymetry data for crustal structure: A prescription for downward continuation | |
| AU2016247039B2 (en) | Geophysical inversion using sparse modeling | |
| Chen et al. | Potential field data interpolation by Taylor series expansion |
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 |