CN104316972A - Strength inversion imaging method for magnetic source gravity - Google Patents

Abstract

The invention relates to a strength inversion imaging method for magnetic source gravity, and the method comprises the steps as follows: changing the magnetic anomaly into the space domain magnetic source gravity anomaly; selecting different cutting radius, adopting interpolation cutting method for separating magnetic source gravity anomaly on magnetic source gravity anomaly of the space domain on different depth space layer, obtaining remaining magnetic source gravity anomaly of each space layer; using iteration method for extending the remaining magnetic source gravity anomaly of each space layer to the top interface of the magnetic body of the space layer, obtaining the magnetic source gravity anomaly of the top interface of the magnetic body; having inversion on the video strength parameter of magnetic source gravity anomaly of the top interface of the magnetic body, obtaining video strength parameter of space layer. The inversion accuracy is high and the inversion multiplicity of solution is reduced, the strength inversion imaging method for magnetic source gravity is widely used for deep mineral product exploration and deep geologic structure study.

Description

A kind of pseudogravity apparent intensity inversion imaging method

Technical field

The present invention relates to a kind of inversion imaging method, particularly about a kind of pseudogravity apparent intensity inversion imaging method in applied geophysics field.

Background technology

At present, two classes are divided into for the method for magnetic field physical property (magnetization) inverting both at home and abroad: linear inversion and non-linear inversion.Concrete algorithm is roughly divided into generalized least square method (Gauss method), method of steepest descent, damped least square method, singular value decomposition method, artificial neural network method BP algorithm, genetic algorithm, simulated annealing and the algorithm etc. based on Radon conversion, multi-scale wavelet edge.Above-mentioned inversion method all achieves certain achievement in theory innovation and the application of small size Real data inversion, but does not perform to due effect in the middle of the actual production for large-area three-dimensional inversion for physical properties.It mainly has problems: 1) when physical parameter and the inverting of the geometric parameter overall situation, need to separate large-scale linear algebraic equation systems, computing time is long; 2) solving equation group mathematics morbid state is very serious, and stability of solution is very poor; 3) the mutual synergistic effect of the positive negative anomaly of magnetic source dipole causes abnormal inversion solution low precision, resolution low.

Summary of the invention

For the problems referred to above, the object of this invention is to provide a kind of inversion accuracy high and the pseudogravity apparent intensity inversion imaging method of inverting multi-solution can be reduced.

For achieving the above object, the present invention takes following technical scheme: a kind of pseudogravity apparent intensity inversion imaging method, it comprises the following steps: gravity anomaly due to magnetic body △ g (x, y) 1) magnetic anomaly △ T (x, y) being converted to spatial domain; 2) choose different cut radius, the gravity anomaly due to magnetic body adopting gravity anomaly due to magnetic body △ g (x, y) of interpolating cut method to spatial domain to carry out different depth space layer is separated, and obtains the residue gravity anomaly due to magnetic body △ g of each space layer _j(x, y, 0); 3) adopt process of iteration by the residue gravity anomaly due to magnetic body △ g of each space layer _j(x, y, 0) downward continuation, on the interface, top of the magnetic of each space layer, obtains gravity anomaly due to magnetic body △ g (x, y, the h at interface, magnetic top _j); 4) to gravity anomaly due to magnetic body △ g (x, y, the h at interface, magnetic top _j) apparent intensity parameter carry out inverting, obtain the apparent intensity parameter of each space layer.

Described step 1) in, magnetic anomaly △ T (x, y) is converted to gravity anomaly due to magnetic body △ g (x, y) of spatial domain, it comprises the following steps: the magnetic anomaly △ T (x, y) that ground observation arrives by (I) converts the magnetic anomaly of frequency field to namely

$\mathrm{\Δ}\tilde{T}(u,v)=F\left[\mathrm{\ΔT}(x,y)\right],$

In formula, u and v represents the wave number in x direction and y direction respectively; (II) by the magnetic anomaly of frequency field converting to pole magnetic anomaly namely

${\tilde{Z}}_{\⊥}(u,v)=\frac{u^2{+v}^2}{[i(P_{0}u+Q_{0}v)+R_0\sqrt{u^2+v^2}]}\·\frac{1}{[i(P_{1}u+Q_{1}v)+R_1\sqrt{u^2+v^2}]}\·\tilde{T}(u,v),$

In formula, i represents imaginary component, (P ₀, Q ₀, R ₀) represent the direction cosine on stratum, magnetic place, (P ₁, Q ₁, R ₁) represent the direction cosine of magnetic on earth's surface, (P ₂, Q ₂, R ₂) represent that magnetic is the direction cosine under the direction of magnetization of 90 ° at magnetic dip; (III) pole magnetic anomaly will be changed dipole magnetic source do class gravity one pole source process, namely right along magnetic south poles direction integral, obtain the gravity anomaly due to magnetic body of frequency field namely

$\mathrm{\Δ}\tilde{g}(u,v)=\frac{\mathrm{G\ρ}}{2\mathrm{\πM}}\frac{\sqrt{u^2+v^2}}{i(P_{0}u+Q_{0}v)+R_0\sqrt{u^2+v^3}}\frac{1}{i(P_{1}u+Q_{1}v)+R_1\sqrt{u^2+v^2}}\·\mathrm{\Δ}\tilde{T}(u,v),$

In formula, M represents magnetization intensity vector, and G represents universal gravitational constant, and ρ represents the residual density of magnetic; (IV) to the gravity anomaly due to magnetic body of frequency field carry out inversefouriertransform, obtain gravity anomaly due to magnetic body △ g (x, y) of spatial domain, namely

$\mathrm{\Δg}(x,y)=F^{-1}\left[\mathrm{\Δ}\tilde{g}(u,v)\right].$

Described step 2) in, obtain the residue gravity anomaly due to magnetic body △ g of each space layer _j(x, y, 0), it specifically comprises the following steps: (I) chooses different cut radius r ₁, r ₂, r ₃r _n, wherein, n represents the number of cut radius, and adopt interpolation cutting method to be separated the gravity anomaly due to magnetic body that gravity anomaly due to magnetic body △ g (x, y) carries out different depth space layer, the shallow source obtained under different cut radius is abnormal:

L ₁(x,y),L ₂(x,y),L ₃(x,y)…L _n(x,y)，

(II) abnormal according to the shallow source under different cut radius, obtain the residue gravity anomaly due to magnetic body of each space layer:

$\left\{\begin{array}{c}\mathrm{\Δ}{g}_0(x,y,0)=L_1(x,y)\\ \mathrm{\Δ}{g}_1(x,y,0)=L_2(x,y)-L_1(x,y)\\ \mathrm{\Δ}{g}_2(x,y,0)=L_3(x,y)-L_2(x,y)\\ .\\ .\\ .\\ \mathrm{\Δ}{g}_j(x,y,0)=L_{j+1}(x,y)-L_j(x,y)\\ .\\ .\\ .\\ \mathrm{\Δ}{g}_{n-1}(x,y,0)=L_n(x,y)-L_{n-1}(x,y)\\ \mathrm{\Δ}{g}_n(x,y,0)=R_n(x,y)\end{array}\right.,$

In formula, △ g _jthe gravity anomaly due to magnetic body that (x, y, 0) expression jth layer magnetic reflects on ground, namely remains gravity anomaly due to magnetic body.

Described step 3) in, adopt process of iteration downward continuation, obtain gravity anomaly due to magnetic body △ g (x, y, the h at interface, magnetic top _j), it specifically comprises the following steps: (I) is by the exceptional value G of s point in the exceptional value △ g (x, y, 0) on the Γ A of ground _asbe put on the vertical projection point s of plane Γ B, as the exceptional value on plane Γ B namely and the frequency spectrum of exceptional value on Γ B is obtained by Fourier transform (II) when not having field source between ground Γ A and plane Γ B, gravity anomaly meets Laplace's equation, adopts upward continuation formula:

$\mathrm{\Δg}(x,y,h)=F^{-1}\left[e^{-2\mathrm{\π}\left|h\right|\sqrt{u^2+v^2}}\mathrm{\Δ}\tilde{g}(u,v,0)\right],$

By the exceptional value on plane Γ B frequency spectrum substitute into the exceptional value that above formula obtains s point on the Γ A of ground that is:

$G_{\mathrm{As}}^{\left(1\right)}(x,y,h)=F^{-1}\left[e^{2\mathrm{\π}\left|h\right|\sqrt{u^2+v^2}}{G}_{\mathrm{Bs}}^{\left(1\right)}(u,v,0)\right],$

In formula, represent the Fourier transform of the gravity anomaly due to magnetic body value △ g (x, y, 0) on ground Γ A (z=0); F ^-1represent inversefouriertransform, h represents continuation height; (III) G is utilized _aswith difference to the exceptional value on plane Γ B revise, obtain new exceptional value that is:

$G_{\mathrm{Bs}}^{\left(2\right)}=G_{\mathrm{Bs}}^{\left(1\right)}+l*(G_{\mathrm{As}}-G_{\mathrm{As}}^{\left(1\right)}),$

In formula, l represents step-length; (IV) repeat step (II) and step (III), obtain following iterative formula:

$G_{\mathrm{Bs}}^{(n+1)}=G_{\mathrm{Bs}}^{\left(n\right)}+l*(G_{\mathrm{As}}-G_{\mathrm{As}}^{\left(n\right)}),$

When $|G_{\mathrm{As}}-G_{\mathrm{As}}^{\left(n\right)}|<\mathrm{\&epsiv;}$ Time, from formula (16) $|G_{\mathrm{Bs}}^{(n+1)}-G_{\mathrm{Bs}}^{\left(n\right)}|<\mathrm{k\&epsiv;},$ ε represent one given level off to zero number, then iterations n gets 20 ~ 50 times; (V) adopt and the method that step (I) ~ (IV) is identical, calculate the gravity anomaly due to magnetic body value of all points on plane Γ B and utilize upward continuation formula:

$\mathrm{\&Delta;g}(x,y,h)=F^{-1}\left[e^{-2\mathrm{\&pi;}\left|h\right|\sqrt{u^2+v^2}}\mathrm{\&Delta;}\tilde{g}(u,v,0)\right],$

△ g (x, y, 0) downward continuation h is obtained as inversefouriertransform _jgravity anomaly due to magnetic body △ g (x, y, the h of the degree of depth _j).

Described step 4) in, obtain the apparent intensity parameter of each space layer, it specifically comprises the following steps: (I) will cause gravity anomaly due to magnetic body △ g (x, y, h _j) irregular magnetic be divided into the different elongated vertical prism of M*N size rule, centroid position; If the center of mass point coordinate of prism is (x _0k, y _0k), its residual density is ρ _k, the bottom surface length of side of prism is respectively a and b, and prism height is dh, then the gravity anomaly △ g that each prism produces at interface, top (x, y) point in frequency field _kthe frequency spectrum of (x, y) for:

$\mathrm{\&Delta;}{\tilde{g}}_k(u,v)=2\mathrm{\&pi;G}{\mathrm{\&rho;}}_k\&CenterDot;\frac{4}{\mathrm{uv}}\sin \left(\frac{\mathrm{au}}{2}\right)\sin \left(\frac{\mathrm{bv}}{2}\right)\&CenterDot;\frac{1}{r}\&CenterDot;(1-e^{-\mathrm{dhr}})\&CenterDot;e^{-i({\mathrm{ux}}_{0k}+{\mathrm{vy}}_{0k})},$

In formula, e represents natural constant, and k represents prism sequence number, k=1,2,3 ..., M*N; (II) frequency spectrum of each prism in the gravity anomaly of interface, top (x, y) some generation is superposed, obtain entire depth layer puts the gravity anomaly △ g (x, y) produced frequency spectrum on interface, top (x, y) for:

$\mathrm{\&Delta;}\tilde{g}(u,v)=2\mathrm{\&pi;G}\&CenterDot;\frac{4}{\mathrm{uv}}\sin \left(\frac{\mathrm{au}}{2}\right)\sin \left(\frac{\mathrm{bv}}{2}\right)\&CenterDot;\frac{1}{r}\&CenterDot;(1-e^{-\mathrm{dhr}}){\mathrm{\&Sigma;}}_{k=1}^{\mathrm{MN}}{\mathrm{\&rho;}}_k\&CenterDot;e^{-i({\mathrm{ux}}_{0k}+{\mathrm{vy}}_k)},$

Order

$\widetilde{\mathrm{\&rho;}}(u,v)={\mathrm{\&Sigma;}}_{k=1}^{\mathrm{MN}}{\mathrm{\&rho;}}_j\&CenterDot;e^{-i({\mathrm{ux}}_{0k}+{\mathrm{vy}}_k)},$

Then gravity anomaly due to magnetic body in the spectral response at interface, top is:

$\mathrm{\&Delta;}\tilde{g}(u,v)=2\mathrm{\&pi;G}\&CenterDot;\frac{4}{\mathrm{uv}}\sin \left(\frac{\mathrm{au}}{2}\right)\sin \left(\frac{\mathrm{bv}}{2}\right)\&CenterDot;\frac{1}{r}\&CenterDot;(1-e^{-\mathrm{dhr}})\&CenterDot;\widetilde{\mathrm{\&rho;}}(u,v),$

Then in frequency domain, the pseudodensity at the interface, residue gravity anomaly due to magnetic body top of space layer is:

$\widetilde{\mathrm{\&rho;}}(u,v)=\frac{\mathrm{\&Delta;}\tilde{g}(u,v)}{2\mathrm{\&pi;G}\&CenterDot;\frac{4}{\mathrm{uv}}\sin \left(\frac{\mathrm{au}}{2}\right)\sin \left(\frac{\mathrm{bv}}{2}\right)\&CenterDot;\frac{1}{r}\&CenterDot;(1-e^{-\mathrm{dhr}})};$

(III) to the pseudodensity at the interface, residue gravity anomaly due to magnetic body top of space layer in frequency domain carry out inversefouriertransform, the residual density value obtaining this depth layer end face each point (x, y) is:

$\mathrm{\&rho;}(x,y)=F^{-1}\left[\frac{\mathrm{\&Delta;}\tilde{g}(u,v)}{2\mathrm{\&pi;G}\&CenterDot;\frac{4}{\mathrm{uv}}\sin \left(\frac{\mathrm{au}}{2}\right)\sin \left(\frac{\mathrm{bv}}{2}\right)\&CenterDot;\frac{1}{r}\&CenterDot;(1-e^{-\mathrm{dhr}})}\right];$

(IV) residual density value ρ (x, y) is substituted into formula

$\mathrm{\&kappa;}(x,y)=\mathrm{\&rho;}(x,y)\&CenterDot;\frac{{\mathrm{\&kappa;}}^{\&prime;}}{{\mathrm{\&rho;}}^{\&prime;}},$

Calculate the apparent magnetic susceptibility κ (x, y) of magnetic anomaly; The apparent magnetic susceptibility κ (x, y) of the magnetic anomaly calculated is substituted into formula

M＝κ(x,y)*T ₀，

Calculate the Windows control center M of magnetic anomaly; In formula, the apparent magnetic susceptibility that κ (x, y) is magnetic anomaly, M is the Windows control center of magnetic anomaly, the pseudodensity that ρ (x, y) is pseudogravity, T ₀for normal magnetic field intensity, the known magnetic susceptibility value needed when κ ' and ρ ' is Poisson formula conversion respectively and density value.

The present invention is owing to taking above technical scheme, and it has the following advantages: 1, existing inversion method is the optimized algorithm based on spatial domain, normally physical parameter and the inverting of the geometric parameter overall situation, and refutation process is complicated, and computational accuracy is low; Compared with existing inversion method, inversion method of the present invention only carries out inverting to physical property intensive parameter, substantially increases inversion accuracy, reduces the multi-solution of inverting.2, existing inversion method demand solution large linear systems, calculator memory consumes huge, and inverse time is long, therefore can not carry out large area three-dimensional imaging inverting fast; Inversion method of the present invention carries out in frequency field, fast operation, and without the need to solving large linear systems, calculator memory consumes little, therefore can carry out on a large scale, three-dimensional Windows control center inversion imaging fast.3, when dipole magnetic source adopts traditional class gravitational inversion method (namely using gravity anomaly inversion method direct inversion magnetic anomaly) in prior art, two different positive negative magnetization parameters can be finally inversed by; And pseudogravity apparent intensity inversion imaging method of the present invention is directly perceived, and the center of magnetic can be reflected.4, existing inversion method is mostly for single model, and for multi-source model inverting, then because magnetic anomaly superposes mutually, abnormal equivalence makes inverting more difficult, and inversion accuracy is poor; And pseudogravity apparent intensity inversion imaging method of the present invention adopts multi-source dipole magnetic source, the one pole superimposed anomaly that can carry out similar gravity anomaly is separated, and the high resolving power Windows control center inverting can carried out based on any deep space magnetic anomaly acquisition of information of magnetic source outer (containing interface, magnetic source top) solves.5, inversion method of the present invention can set up three-dimensional Windows control center property distribution.Based on above advantage, the present invention can be widely used in deep ore prospecting and Deep Geological Structures research.

Accompanying drawing explanation

Fig. 1 is the process flow diagram of pseudogravity apparent intensity inversion imaging method of the present invention;

Fig. 2 is the schematic diagram adopting process of iteration downward continuation;

Fig. 3 is theoretical model space distribution and floor projection figure; Wherein, figure (a) is theoretical model spatial distribution map, and figure (b) is the floor projection figure of theoretical model;

Fig. 4 is theoretical model ground magnetic anomaly schematic diagram;

Fig. 5 is that ground magnetic anomaly normalizing is to the abnormal schematic diagram of terrestrial magnetic pole;

Fig. 6 is Model A terrestrial gravity anomaly due to magnetic body schematic diagram;

Fig. 7 is 0 ~ 200 meter of space layer residue magnetic source gravity anomaly chart;

Fig. 8 is 200 ~ 500 meters of space layer residue magnetic source gravity anomaly chart;

Fig. 9 is 200 meters of Shen Ding interface gravity anomaly due to magnetic body figure;

Figure 10 is 500 meters of Shen Ding interface gravity anomaly due to magnetic body figure;

Figure 11 is 200 meters of Shen Ding interface pseudogravity pseudodensity inversion charts;

Figure 12 is 500 meters of Shen Ding interface pseudogravity pseudodensity inversion charts;

Figure 13 is 200 meters of Shen Ding interface pseudogravity Windows control center inversion charts;

Figure 14 is 500 meters of Shen Ding interface pseudogravity Windows control center inversion charts.

Embodiment

Below in conjunction with drawings and Examples, the present invention is described in detail.

As shown in Figure 1, pseudogravity apparent intensity inversion imaging method of the present invention comprises the following steps:

1) magnetic anomaly △ T (x, y) is converted to gravity anomaly due to magnetic body △ g (x, y) of spatial domain;

When there is density difference XOR magnetic contrast between underground geologic bodies, ground can observe its corresponding gravity anomaly △ g (x, or magnetic anomaly △ T (x y), y), because heavy magnetic field meets Laplce's potential field theory, therefore mathematically can do mutual conversion, namely gravity anomaly can be converted into gravity source magnetic anomaly, magnetic anomaly by the form of gradient and also can become gravity anomaly due to magnetic body with the formal transformation of integration.Illustrate the implementation procedure that magnetic anomaly converts gravity anomaly due to magnetic body to below:

(I) the magnetic anomalies △ T (x, y) arrived by ground observation converts the magnetic anomaly of frequency field to

Fourier transform is carried out to the magnetic anomaly △ T (x, y) that ground observes, obtains the frequency spectrum on ground magnetic anomaly △ T (x, y) namely

$\mathrm{\&Delta;}\tilde{T}(u,v)=F\left[\mathrm{\&Delta;T}(x,y)\right]---\left(1\right)$

In formula, u and v represents the wave number in x direction and y direction respectively.

(II) by the magnetic anomaly of frequency field be converted into the magnetic anomaly with gravity direction, i.e. converting to pole magnetic anomaly

Actual measurement magnetic anomalies in setting earth's surface is its initial magnetic dip of magnetic of △ T (x, y) is I ₁, magnetic declination is D ₁, the direction cosine of this magnetic on earth's surface are (P ₁, Q ₁, R ₁); The earth's magnetic dip angle on this stratum, magnetic place is I ₀, geomagnetic declination is D ₀, the direction cosine on this stratum, magnetic place are (P ₀, Q ₀, R ₀); The direction of magnetization of this magnetic is I from magnetic dip ₁being transformed into magnetic dip is I ₂time, magnetic declination is D ₂, the direction cosine of this magnetic under this direction of magnetization are (P ₂, Q ₂, R ₂), be then I by direction of magnetization from magnetic dip ₁converting magnetic dip to is I ₂time, direction of magnetization conversion formula is:

${\tilde{Z}}_a(u,v)=\frac{\sqrt{u^2{+v}^2}}{[i(P_{0}u+Q_{0}v)+R_0\sqrt{u^2+v^2}]}\&CenterDot;\frac{[i(P_{2}u+Q_{2}v)+R_2\sqrt{u^2+v^2}]}{[i(P_{1}u+Q_{1}v)+R_1\sqrt{u^2+v^2}]}\&CenterDot;\tilde{T}(u,v)---\left(2\right)$

In formula, i represents imaginary component, i.e. i ²=-1.This magnetic is at the direction cosine (P on earth's surface ₁, Q ₁, R ₁) be:

$\left\{\begin{array}{c}{P}_1=\cos I_1\cos D_1\\ Q_1=\cos I_1\sin D_1\\ R_1=\sin I_1\end{array}\right.,$

Direction cosine (the P on this stratum, magnetic place ₀, Q ₀, R ₀) be:

$\left\{\begin{array}{c}{P}_0=\cos I_0\cos D_0\\ Q_0=\cos I_0\sin D_0\\ R_0=\sin I_0\end{array}\right.,$

This magnetic is I at magnetic dip ₂direction of magnetization under direction cosine (P ₂, Q ₂, R ₂) be:

$\left\{\begin{array}{c}{P}_2=\cos I_2\cos D_2\\ Q_2=\cos I_2\sin D_2\\ R_2=\sin I_2\end{array}\right.---\left(3\right)$

When the direction of magnetization of this magnetic is I from magnetic dip ₁when being transformed into perpendicular magnetization direction (changing to terrestrial magnetic pole), i.e. I ₂=90 °, by I ₂=90 ° substitute in formula (3), then this magnetic is I at magnetic dip ₂direction of magnetization under direction cosine (P ₂, Q ₂, R ₂) be:

$\left\{\begin{array}{c}{P}_2=0\\ Q_2=0\\ R_2=1\end{array}\right.---\left(4\right)$

Formula (4) brought in formula (2), obtaining is to the frequency spectrum of terrestrial magnetic pole for:

${\tilde{Z}}_{\&perp;}(u,v)=\frac{u^2{+v}^2}{[i(P_{0}u+Q_{0}v)+R_0\sqrt{u^2+v^2}]}\&CenterDot;\frac{1}{[i(P_{1}u+Q_{1}v)+R_1\sqrt{u^2+v^2}]}\&CenterDot;\tilde{T}(u,v)---\left(5\right)$

(III) pole magnetic anomaly will be changed dipole magnetic source do class gravity one pole source process, namely right along magnetic south poles direction integral, obtain the gravity anomaly due to magnetic body of frequency field

Due to magnetic field intensity and the gravity field strength value of the corresponding gradient direction of heavy magnetic potential gradient value difference, namely

$\left\{\begin{array}{c}\frac{\&PartialD;U}{\&PartialD;z}={\tilde{Z}}_{\&perp;}(u,v)\\ \frac{\&PartialD;V}{\&PartialD;z}=\mathrm{\&Delta;}\tilde{g}(u,v)\end{array}\right.---\left(6\right)$

And the relation between magnetic potential U and gravitation position V meets Poisson formula, namely

$U=-\frac{1}{4\mathrm{\&pi;G\&rho;}}M\&CenterDot;\mathrm{grad}V---\left(7\right)$

In formula, M represents magnetization intensity vector, and G represents universal gravitational constant, and ρ represents the residual density of magnetic.

The gravity anomaly due to magnetic body of frequency field is calculated by formula (5) ~ formula (7) for:

$\mathrm{\&Delta;}\tilde{g}(u,v)=\frac{\mathrm{G\&rho;}}{2\mathrm{\&pi;M}}\frac{\sqrt{u^2+v^2}}{i(P_{0}u+Q_{0}v)+R_0\sqrt{u^2+v^3}}\frac{1}{i(P_{1}u+Q_{1}v)+R_1\sqrt{u^2+v^2}}\&CenterDot;\mathrm{\&Delta;}\tilde{T}(u,v)---\left(8\right)$

(IV) to the gravity anomaly due to magnetic body of frequency field carry out inversefouriertransform, obtain gravity anomaly due to magnetic body △ g (x, y) of spatial domain, namely

$\mathrm{\&Delta;g}(x,y)=F^{-1}\left[\mathrm{\&Delta;}\tilde{g}(u,v)\right]---\left(9\right)$

2) choose different cut radius, the gravity anomaly due to magnetic body adopting gravity anomaly due to magnetic body △ g (x, y) of interpolating cut method to spatial domain to carry out different depth space layer is separated, and obtains the residue gravity anomaly due to magnetic body △ g of each space layer _j(x, y, 0), it specifically comprises the following steps:

(I) different cut radius r is chosen ₁, r ₂, r ₃r _n, wherein, n represents the number of cut radius, and adopt interpolation cutting method to be separated the gravity anomaly due to magnetic body that gravity anomaly due to magnetic body △ g (x, y) carries out different depth space layer, the shallow source obtained under different cut radius is abnormal:

L ₁(x,y),L ₂(x,y),L ₃(x,y)…L _n(x,y) (10)

(II) abnormal according to the shallow source under different cut radius, obtain the residue gravity anomaly due to magnetic body of each space layer:

$\left\{\begin{array}{c}\mathrm{\&Delta;}{g}_0(x,y,0)=L_1(x,y)\\ \mathrm{\&Delta;}{g}_1(x,y,0)=L_2(x,y)-L_1(x,y)\\ \mathrm{\&Delta;}{g}_2(x,y,0)=L_3(x,y)-L_2(x,y)\\ .\\ .\\ .\\ \mathrm{\&Delta;}{g}_j(x,y,0)=L_{j+1}(x,y)-L_j(x,y)\\ .\\ .\\ .\\ \mathrm{\&Delta;}{g}_{n-1}(x,y,0)=L_n(x,y)-L_{n-1}(x,y)\\ \mathrm{\&Delta;}{g}_n(x,y,0)=R_n(x,y)\end{array}\right.---\left(11\right)$

In formula, △ g _j(x, y, 0) represents jth (j=0,1,2,3 ... n) gravity anomaly due to magnetic body that reflects on ground of layer magnetic, namely remains gravity anomaly due to magnetic body.

3) adopt process of iteration by the residue gravity anomaly due to magnetic body △ g of each space layer _j(x, y, 0) downward continuation, on the interface, top of the magnetic of each space layer, obtains gravity anomaly due to magnetic body △ g (x, y, the h at interface, magnetic top _j), wherein, △ g (x, y, h _j) represent j layer residue gravity anomaly due to magnetic body △ g _j(x, y, 0) downward continuation is to below ground h _jthe degree of depth, obtains the gravity anomaly due to magnetic body at interface, magnetic top, degree of depth h _jfor:

$\left\{\begin{array}{c}{h}_0=0,\\ h_1=r_1,\\ .\\ .\\ .,\\ h_j=r_j,\\ .\\ .\\ .,\\ h_n=r_n\end{array}\right.---\left(12\right)$

Adopt process of iteration downward continuation, obtain gravity anomaly due to magnetic body △ g (x, y, the h at interface, magnetic top _j) concrete steps comprise:

(I) as shown in Figure 2, by the exceptional value G of s point in the exceptional value △ g (x, y, 0) on the Γ A of ground _asbe put on the vertical projection point s of plane Γ B, as the exceptional value on plane Γ B namely and the frequency spectrum of exceptional value on Γ B is obtained by Fourier transform

(II) when not having field source between ground Γ A and plane Γ B, gravity anomaly meets Laplace's equation, adopts upward continuation formula:

$\mathrm{\&Delta;g}(x,y,h)=F^{-1}\left[e^{-2\mathrm{\&pi;}\left|h\right|\sqrt{u^2+v^2}}\mathrm{\&Delta;}\tilde{g}(u,v,0)\right]---\left(13\right)$

By the exceptional value on plane Γ B frequency spectrum substitute into the exceptional value that above formula obtains s point on the Γ A of ground that is:

$G_{\mathrm{As}}^{\left(1\right)}(x,y,h)=F^{-1}\left[e^{2\mathrm{\&pi;}\left|h\right|\sqrt{u^2+v^2}}{G}_{\mathrm{Bs}}^{\left(1\right)}(u,v,0)\right]---\left(14\right)$

In formula (13), represent the Fourier transform of the gravity anomaly due to magnetic body value △ g (x, y, 0) on ground Γ A (z=0); F ^-1represent inversefouriertransform, h represents continuation height.

(III) G is utilized _aswith difference to the exceptional value on plane Γ B revise, obtain new exceptional value that is:

$G_{\mathrm{Bs}}^{\left(2\right)}=G_{\mathrm{Bs}}^{\left(1\right)}+l*(G_{\mathrm{As}}-G_{\mathrm{As}}^{\left(1\right)})---\left(15\right)$

In formula, l represents step-length, and general value is 1.

(IV) repeat step (II) and step (III), obtain following iterative formula:

$G_{\mathrm{Bs}}^{(n+1)}=G_{\mathrm{Bs}}^{\left(n\right)}+l*(G_{\mathrm{As}}-G_{\mathrm{As}}^{\left(n\right)})---\left(16\right)$

When $|G_{\mathrm{As}}-G_{\mathrm{As}}^{\left(n\right)}|<\mathrm{\&epsiv;}$ Time, from formula (16) $|G_{\mathrm{Bs}}^{(n+1)}-G_{\mathrm{Bs}}^{\left(n\right)}|<\mathrm{k\&epsiv;},$ ε represent one given level off to zero number, then $G_{\mathrm{Bs}}^{(n+1)}\&ap;G_{\mathrm{Bs}}^{\left(n\right)}.$

General iterations n gets 20 ~ 50 times.

(V) adopt and the method that step (I) ~ (IV) is identical, calculate the gravity anomaly due to magnetic body value of all points on plane Γ B and utilize upward continuation formula (13) to obtain gravity anomaly due to magnetic body △ g (x, y, the h of △ g (x, y, 0) the downward continuation hj degree of depth _j).

4) to interface, top gravity anomaly due to magnetic body △ g (x, y, h _j) apparent intensity parameter carry out inverting, obtain the apparent intensity parameter of each space layer, it specifically comprises the following steps:

(I) gravity anomaly due to magnetic body △ g (x, y, h will be caused _j) irregular magnetic be divided into the different elongated vertical prism of M*N size rule, centroid position.

If the center of mass point coordinate of prism is (x _0k, y _0k), its residual density is ρ _k(k=1,2,3 ..., M*N), the bottom surface length of side of prism is respectively a and b, and prism height is dh, then in frequency field, the frequency spectrum of the gravity anomaly △ gk (x, y) produced put by each prism on interface, top (x, y) for:

$\mathrm{\&Delta;}{\tilde{g}}_k(u,v)=2\mathrm{\&pi;G}{\mathrm{\&rho;}}_k\&CenterDot;\frac{4}{\mathrm{uv}}\sin \left(\frac{\mathrm{au}}{2}\right)\sin \left(\frac{\mathrm{bv}}{2}\right)\&CenterDot;\frac{1}{r}\&CenterDot;(1-e^{-\mathrm{dhr}})\&CenterDot;e^{-i({\mathrm{ux}}_{0k}+{\mathrm{vy}}_{0k})}---\left(17\right)$

In formula, e represents natural constant, and k represents prism sequence number.

(II) frequency spectrum of each prism in the gravity anomaly of interface, top (x, y) some generation is superposed, obtain entire depth layer puts the gravity anomaly △ g (x, y) produced frequency spectrum on interface, top (x, y) for:

$\mathrm{\&Delta;}\tilde{g}(u,v)=2\mathrm{\&pi;G}\&CenterDot;\frac{4}{\mathrm{uv}}\sin \left(\frac{\mathrm{au}}{2}\right)\sin \left(\frac{\mathrm{bv}}{2}\right)\&CenterDot;\frac{1}{r}\&CenterDot;(1-e^{-\mathrm{dhr}}){\mathrm{\&Sigma;}}_{k=1}^{\mathrm{MN}}{\mathrm{\&rho;}}_k\&CenterDot;e^{-i({\mathrm{ux}}_{0k}+{\mathrm{vy}}_k)}---\left(18\right)$

Order

$\widetilde{\mathrm{\&rho;}}(u,v)={\mathrm{\&Sigma;}}_{k=1}^{\mathrm{MN}}{\mathrm{\&rho;}}_j\&CenterDot;e^{-i({\mathrm{ux}}_{0k}+{\mathrm{vy}}_k)}---\left(19\right)$

Then gravity anomaly due to magnetic body in the spectral response at interface, top is:

$\mathrm{\&Delta;}\tilde{g}(u,v)=2\mathrm{\&pi;G}\&CenterDot;\frac{4}{\mathrm{uv}}\sin \left(\frac{\mathrm{au}}{2}\right)\sin \left(\frac{\mathrm{bv}}{2}\right)\&CenterDot;\frac{1}{r}\&CenterDot;(1-e^{-\mathrm{dhr}})\&CenterDot;\widetilde{\mathrm{\&rho;}}(u,v)---\left(20\right)$

The pseudodensity being obtained interface, gravity anomaly top, space layer remanence source in frequency domain by formula (20) is:

$\widetilde{\mathrm{\&rho;}}(u,v)=\frac{\mathrm{\&Delta;}\tilde{g}(u,v)}{2\mathrm{\&pi;G}\&CenterDot;\frac{4}{\mathrm{uv}}\sin \left(\frac{\mathrm{au}}{2}\right)\sin \left(\frac{\mathrm{bv}}{2}\right)\&CenterDot;\frac{1}{r}\&CenterDot;(1-e^{-\mathrm{dhr}})}---\left(21\right)$

(III) to the pseudodensity at interface, gravity anomaly top, space layer remanence source in frequency domain carry out inversefouriertransform, the residual density value obtaining this depth layer end face each point (x, y) is:

$\mathrm{\&rho;}(x,y)=F^{-1}\left[\frac{\mathrm{\&Delta;}\tilde{g}(u,v)}{2\mathrm{\&pi;G}\&CenterDot;\frac{4}{\mathrm{uv}}\sin \left(\frac{\mathrm{au}}{2}\right)\sin \left(\frac{\mathrm{bv}}{2}\right)\&CenterDot;\frac{1}{r}\&CenterDot;(1-e^{-\mathrm{dhr}})}\right]---\left(22\right)$

(IV) according to Poisson formula, there is following linear relationship in magnetic susceptibility and density:

$\mathrm{\&kappa;}(x,y)=\mathrm{\&rho;}(x,y)\&CenterDot;\frac{{\mathrm{\&kappa;}}^{\&prime;}}{{\mathrm{\&rho;}}^{\&prime;}}---\left(23\right)$

M＝κ(x,y)*T ₀ (24)

In formula, the apparent magnetic susceptibility that κ (x, y) is magnetic anomaly, M is the Windows control center of magnetic anomaly, the pseudodensity that ρ (x, y) is pseudogravity, T ₀for normal magnetic field intensity, the known magnetic susceptibility value needed when κ ' and ρ ' is Poisson formula conversion respectively and density value.

By the residual density value ρ (x, y) calculated, the i.e. pseudodensity of pseudogravity, substitute in formula (23), calculate the apparent magnetic susceptibility κ (x, y) of magnetic anomaly.The apparent magnetic susceptibility κ (x, y) of the magnetic anomaly calculated is substituted in formula (24), calculates the Windows control center M of magnetic anomaly.

Embodiment 1: in order to describe the process of inversion imaging method of the present invention and technical advantage in detail, creates a type cylinder models group as shown in Figure 3, each model post just drill geometric parameter and physical parameter is as shown in table 1.

The design parameter of table 1 model

Adopt inversion imaging method of the present invention to carry out inverting to the geologic body of this model representation, its specific implementation process is:

(1) obtain ground magnetic field value as shown in Figure 4, and regard it as ground observation data;

(2) gravity direction conversion is done to ground observation magnetic anomaly data, obtain change pole magnetic anomaly as shown in Figure 5;

(3) conversion of one pole source is done to change pole dipole magnetic source, obtain gravity anomaly due to magnetic body as shown in Figure 6;

(4) utilize the cut radius of 200 meters and 500 meters to carry out inter-subspecific hybrid to ground pseudogravity respectively, obtain the ground residue gravity anomaly due to magnetic body of 0 ~ 200 meter of space layer as shown in Figure 7 and 200 ~ 500 meters of space layer as shown in Figure 8;

(5) process of iteration is adopted to carry out downward continuation to 0 ~ 200 meter of space layer residue pseudogravity, obtain interface, 200 meters of dark magnetic tops gravity anomaly due to magnetic body as shown in Figure 9, adopt process of iteration to carry out downward continuation to 200 ~ 500 meters of space layer residue pseudogravity, obtain interface, 500 meters of dark magnetic tops gravity anomaly due to magnetic body as shown in Figure 10;

(6) respectively pseudodensity inverting is done to 200 meters of dark magnetic top interface gravity anomaly due to magnetic body and interface, 500 meters of dark magnetic tops pseudogravity, acquisition 200 meters of Shen Ding interface pseudogravity pseudodensity inversion charts as shown in figure 11 and 500 meters of Shen Ding interface pseudogravity pseudodensity inversion charts as shown in figure 12;

(7) respectively Windows control center conversion is done to the pseudodensity of 200 meters of dark pseudodensities and 500 meters of dark magnetics, obtain 200 meters of Shen Ding interface pseudogravity apparent intensity inversion charts and 500 meters of Shen Ding interface pseudogravity apparent intensity inversion charts as shown in figure 14 as shown in fig. 13 that.

As can be seen from Fig. 4 ~ Figure 14, adopt pseudogravity apparent intensity formation method of the present invention, the parameter as shown in Figure 13 and Figure 14 that inverting obtains, compared with model parameter as shown in Figure 3, has good degree of accuracy.Wherein, the deep focus inverting magnetization is 75, with gross data 80 relatively, shallow source magnetization inverting value is 100, and coincide with theoretical value 100, this illustration method theory is feasible.In actual applications, owing to there is certain randomness when choosing cut radius, therefore the interlayer residual anomaly separated truly can not reflect the locus of actual magnetic gonosome, adopt " interface, magnetic top " after process of iteration downward continuation may not overlap with the interface, top of true magnetic, make inverting intensive parameter out and actual parameter value there is certain difference.But inverted parameters does not affect the qualitative comparison of magnetic of laterally multiple magnetic, and therefore, the parameter of institute's inverting is called apparent intensity parametric inversion by the present invention.

The various embodiments described above are only for illustration of the present invention; wherein the structure of each parts, connected mode and method step etc. all can change to some extent; every equivalents of carrying out on the basis of technical solution of the present invention and improvement, all should not get rid of outside protection scope of the present invention.

Claims (5)

1. a pseudogravity apparent intensity inversion imaging method, it comprises the following steps:

1) magnetic anomaly △ T (x, y) is converted to gravity anomaly due to magnetic body △ g (x, y) of spatial domain;

2) choose different cut radius, the gravity anomaly due to magnetic body adopting gravity anomaly due to magnetic body △ g (x, y) of interpolating cut method to spatial domain to carry out different depth space layer is separated, and obtains the residue gravity anomaly due to magnetic body △ g of each space layer _j(x, y, 0);

3) adopt process of iteration by the residue gravity anomaly due to magnetic body △ g of each space layer _j(x, y, 0) downward continuation, on the interface, top of the magnetic of each space layer, obtains gravity anomaly due to magnetic body △ g (x, y, the h at interface, magnetic top _j);

4) to gravity anomaly due to magnetic body △ g (x, y, the h at interface, magnetic top _j) apparent intensity parameter carry out inverting, obtain the apparent intensity parameter of each space layer.

2. a kind of pseudogravity apparent intensity inversion imaging method as claimed in claim 1, it is characterized in that: described step 1) in, magnetic anomaly △ T (x, y) is converted to the gravity anomaly due to magnetic body △ g (x of spatial domain, y), it comprises the following steps:

(I) the magnetic anomaly △ T (x, y) arrived by ground observation converts the magnetic anomaly of frequency field to namely

$\mathrm{\&Delta;}\tilde{T}(u,v)=F\left[\mathrm{\&Delta;T}(x,y)\right],$

In formula, u and v represents the wave number in x direction and y direction respectively;

(II) by the magnetic anomaly of frequency field converting to pole magnetic anomaly namely

${\tilde{Z}}_{\&perp;}(u,v)=\frac{u^2+v^2}{[i(P_{0}u+Q_{0}v)+R_0\sqrt{u^2+v^2}]}\&CenterDot;\frac{1}{[i(P_{1}u+Q_{1}v)+R_1\sqrt{u^2+v^2}]}\&CenterDot;\tilde{T}(u,v),$

In formula, i represents imaginary component, (P ₀, Q ₀, R ₀) represent the direction cosine on stratum, magnetic place, (P ₁, Q ₁, R ₁) represent the direction cosine of magnetic on earth's surface, (P ₂, Q ₂, R ₂) represent that magnetic is the direction cosine under the direction of magnetization of 90 ° at magnetic dip;

(III) pole magnetic anomaly will be changed dipole magnetic source do class gravity one pole source process, namely right along magnetic south poles direction integral, obtain the gravity anomaly due to magnetic body of frequency field namely

$\mathrm{\&Delta;}\tilde{g}(u,v)=\frac{\mathrm{G\&rho;}}{2\mathrm{\&pi;M}}\frac{\sqrt{u^2+v^2}}{i(P_{0}u+Q_{0}v)+R_0\sqrt{u^2+v^2}}\frac{1}{i(P_{1}u+Q_{1}v)+R_1\sqrt{u^2+v^2}}\&CenterDot;\mathrm{\&Delta;}\tilde{T}(u,v),$

In formula, M represents magnetization intensity vector, and G represents universal gravitational constant, and ρ represents the residual density of magnetic;

(IV) to the gravity anomaly due to magnetic body of frequency field carry out inversefouriertransform, obtain gravity anomaly due to magnetic body △ g (x, y) of spatial domain, namely

$\mathrm{\&Delta;g}(x,y)=F^{-1}\left[\mathrm{\&Delta;}\tilde{g}(u,v)\right].$

3. a kind of pseudogravity apparent intensity inversion imaging method as claimed in claim 1 or 2, is characterized in that: described step 2) in, obtain the residue gravity anomaly due to magnetic body △ g of each space layer _j(x, y, 0), it specifically comprises the following steps:

(I) different cut radius r is chosen ₁, r ₂, r ₃r _n, wherein, n represents the number of cut radius, and adopt interpolation cutting method to be separated the gravity anomaly due to magnetic body that gravity anomaly due to magnetic body △ g (x, y) carries out different depth space layer, the shallow source obtained under different cut radius is abnormal:

L ₁(x,y),L ₂(x,y),L ₃(x,y)…L _n(x,y)，

(II) abnormal according to the shallow source under different cut radius, obtain the residue gravity anomaly due to magnetic body of each space layer:

$\left\{\begin{array}{c}{\mathrm{\&Delta;g}}_0(x,y,0)=L_1(x,y)\\ {\mathrm{\&Delta;g}}_1(x,y,0)=L_2(x,y)-L_1(x,y)\\ {\mathrm{\&Delta;g}}_2(x,y,0)=L_3(x,y)-L_2(x,y)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ {\mathrm{\&Delta;g}}_j(x,y,0)=L_{j+1}(x,y)-L_j(x,y)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ {\mathrm{\&Delta;g}}_{n-1}(x,y,0)=L_n(x,y)-L_{n-1}(x,y)\\ {\mathrm{\&Delta;g}}_n(x,y,0)=R_n(x,y)\end{array}\right.,$

In formula, △ g _jthe gravity anomaly due to magnetic body that (x, y, 0) expression jth layer magnetic reflects on ground, namely remains gravity anomaly due to magnetic body.

4. a kind of pseudogravity apparent intensity inversion imaging method as claimed in claim 1 or 2, is characterized in that: described step 3) in, adopt process of iteration downward continuation, obtain gravity anomaly due to magnetic body △ g (x, y, the h at interface, magnetic top _j), it specifically comprises the following steps:

(I) by the exceptional value G of s point in the exceptional value △ g (x, y, 0) on the Γ A of ground _asbe put on the vertical projection point s of plane Γ B, as the exceptional value on plane Γ B namely and the frequency spectrum of exceptional value on Γ B is obtained by Fourier transform

(II) when not having field source between ground Γ A and plane Γ B, gravity anomaly meets Laplace's equation, adopts upward continuation formula:

$\mathrm{\&Delta;g}(x,y,h)=F^{-1}\left[e^{-2\mathrm{\&pi;}\left|h\right|\sqrt{u^2+v^2}}\mathrm{\&Delta;}\tilde{g}(u,v,0)\right],$

By the exceptional value on plane Γ B frequency spectrum substitute into the exceptional value that above formula obtains s point on the Γ A of ground that is:

$G_{\mathrm{As}}^{\left(1\right)}(x,y,h)=F^{-1}\left[e^{-2\mathrm{\&pi;}\left|h\right|\sqrt{u^2+v^2}}{G}_{\mathrm{Bs}}^{\left(1\right)}(u,v,0)\right],$

In formula, represent the Fourier transform of the gravity anomaly due to magnetic body value △ g (x, y, 0) on ground Γ A (z=0); F ^-1represent inversefouriertransform, h represents continuation height;

(III) G is utilized _aswith difference to the exceptional value on plane Γ B revise, obtain new exceptional value that is:

$G_{\mathrm{Bs}}^{\left(2\right)}=G_{\mathrm{Bs}}^{\left(1\right)}+l*(G_{\mathrm{As}}-G_{\mathrm{As}}^{\left(1\right)}),$

In formula, l represents step-length;

(IV) repeat step (II) and step (III), obtain following iterative formula:

$G_{\mathrm{Bs}}^{(n+1)}=G_{\mathrm{Bs}}^{\left(n\right)}+l*(G_{\mathrm{As}}-G_{\mathrm{As}}^{\left(n\right)}),$

When time, from formula (16) ε represent one given level off to zero number, then iterations n gets 20 ~ 50 times;

(V) adopt and the method that step (I) ~ (IV) is identical, calculate the gravity anomaly due to magnetic body value of all points on plane Γ B and utilize upward continuation formula:

$\mathrm{\&Delta;g}(x,y,h)=F^{-1}\left[e^{-2\mathrm{\&pi;}\left|h\right|\sqrt{u^2+v^2}}\mathrm{\&Delta;}\tilde{g}(u,v,0)\right],$

△ g (x, y, 0) downward continuation h is obtained as inversefouriertransform _jgravity anomaly due to magnetic body △ g (x, y, the h of the degree of depth _j).

5. a kind of pseudogravity apparent intensity inversion imaging method as claimed in claim 1 or 2, is characterized in that: described step 4) in, obtain the apparent intensity parameter of each space layer, it specifically comprises the following steps:

(I) gravity anomaly due to magnetic body △ g (x, y, h will be caused _j) irregular magnetic be divided into the different elongated vertical prism of M*N size rule, centroid position;

If the center of mass point coordinate of prism is (x _0k, y _0k), its residual density is ρ _k, the bottom surface length of side of prism is respectively a and b, and prism height is dh, then the gravity anomaly △ g that each prism produces at interface, top (x, y) point in frequency field _kthe frequency spectrum of (x, y) for:

$\mathrm{\&Delta;}{\tilde{g}}_k(u,v)=2\mathrm{\&pi;G}{\mathrm{\&rho;}}_k\&CenterDot;\frac{4}{\mathrm{uv}}\sin \left(\frac{\mathrm{au}}{2}\right)\sin \left(\frac{\mathrm{bv}}{2}\right)\&CenterDot;\frac{1}{r}\&CenterDot;(1-e^{-\mathrm{dhr}})\&CenterDot;e^{-i({\mathrm{ux}}_{0k}+{\mathrm{vy}}_{0k})},$

In formula, e represents natural constant, and k represents prism sequence number, k=1,2,3 ..., M*N;

(II) frequency spectrum of each prism in the gravity anomaly of interface, top (x, y) some generation is superposed, obtain entire depth layer puts the gravity anomaly △ g (x, y) produced frequency spectrum on interface, top (x, y) for:

$\mathrm{\&Delta;}\tilde{g}(u,v)=2\mathrm{\&pi;G}\&CenterDot;\frac{4}{\mathrm{uv}}\sin \left(\frac{\mathrm{au}}{2}\right)\sin \left(\frac{\mathrm{bv}}{2}\right)\&CenterDot;\frac{1}{r}\&CenterDot;(1-e^{-\mathrm{dhr}}){\mathrm{\&Sigma;}}_{k=1}^{\mathrm{MN}}{\mathrm{\&rho;}}_k\&CenterDot;e^{-i({\mathrm{ux}}_{0k}+{\mathrm{vy}}_k)},$

Order

$\widetilde{\mathrm{\&rho;}}(u,v)={\mathrm{\&Sigma;}}_{k=1}^{\mathrm{MN}}{\mathrm{\&rho;}}_j\&CenterDot;e^{-i({\mathrm{ux}}_{0k}+{\mathrm{vy}}_k)},$

Then gravity anomaly due to magnetic body in the spectral response at interface, top is:

$\mathrm{\&Delta;}\tilde{g}(u,v)=2\mathrm{\&pi;G}\&CenterDot;\frac{4}{\mathrm{uv}}\sin \left(\frac{\mathrm{au}}{2}\right)\sin \left(\frac{\mathrm{by}}{2}\right)\&CenterDot;\frac{1}{r}\&CenterDot;(1-e^{-\mathrm{dhr}})\&CenterDot;\widetilde{\mathrm{\&rho;}}(u,v),$

Then in frequency domain, the pseudodensity at the interface, residue gravity anomaly due to magnetic body top of space layer is:

$\widetilde{\mathrm{\&rho;}}(u,v)=\frac{\mathrm{\&Delta;}\tilde{g}(u,v)}{2\mathrm{\&pi;G}\&CenterDot;\frac{4}{\mathrm{uv}}\sin \left(\frac{\mathrm{au}}{2}\right)\sin \left(\frac{\mathrm{by}}{2}\right)\&CenterDot;\frac{1}{r}\&CenterDot;(1-e^{-\mathrm{dhr}})};$

(III) to the pseudodensity at the interface, residue gravity anomaly due to magnetic body top of space layer in frequency domain carry out inversefouriertransform, the residual density value obtaining this depth layer end face each point (x, y) is:

$\mathrm{\&rho;}(x,y)=F^{-1}\left[\frac{\mathrm{\&Delta;}\tilde{g}(u,v)}{2\mathrm{\&pi;G}\&CenterDot;\frac{4}{\mathrm{uv}}\sin \left(\frac{\mathrm{au}}{2}\right)\sin \left(\frac{\mathrm{bv}}{2}\right)\&CenterDot;\frac{1}{r}\&CenterDot;(1-e^{-\mathrm{dhr}})}\right];$

(IV) residual density value ρ (x, y) is substituted into formula

$\mathrm{\&kappa;}(x,y)=\mathrm{\&rho;}(x,y)\&CenterDot;\frac{{\mathrm{\&kappa;}}^{\&prime;}}{{\mathrm{\&rho;}}^{\&prime;}},$

Calculate the apparent magnetic susceptibility κ (x, y) of magnetic anomaly; The apparent magnetic susceptibility κ (x, y) of the magnetic anomaly calculated is substituted into formula

M＝κ(x,y)*T ₀，

Calculate the Windows control center M of magnetic anomaly;

In formula, the apparent magnetic susceptibility that κ (x, y) is magnetic anomaly, M is the Windows control center of magnetic anomaly, the pseudodensity that ρ (x, y) is pseudogravity, T ₀for normal magnetic field intensity, the known magnetic susceptibility value needed when κ ' and ρ ' is Poisson formula conversion respectively and density value.