阅读说明：本技术 一种基于重磁界面反演的大地热流估计方法 (A kind of terrestrial heat flow estimation method based on boundary surface of gravitation and magnetism inverting ) 是由 曾昭发 王卓 刘卓 李静 于 2019-07-03 设计创作，主要内容包括：本发明公开了一种基于重磁界面反演的大地热流估计方法,该方法包括(S1)利用改进后的Parker-Oldenburg界面反演算法对重磁数据进行反演得出居里面和莫霍面,所述莫霍面和居里面的反演算法分别为：<Image he="57" wi="700" file="DDA0002116927700000011.GIF" imgContent="drawing" imgFormat="GIF" orientation="portrait" inline="no"></Image><Image he="54" wi="700" file="DDA0002116927700000012.GIF" imgContent="drawing" imgFormat="GIF" orientation="portrait" inline="no"></Image>其中,h<Sub>0</Sub>为上下界面的平均深度,Δh<Sub>2</Sub>和ΔH<Sub>1</Sub>为上下界面与平均深度的差值,s为迭代步长,ρ<Sub>0</Sub>为地表地质介质的剩余密度,ω为圆波数,F[]为重力异常的傅里叶变换,Δg为重力异常,G为万有引力常数,a为密度随深度变化指数；Δz为磁异常,μ<Sub>0</Sub>为真空中磁导率,b为磁化率随深度变化指数；(S2)利用居里面和莫霍面对地热靶区进行评价。通过上述方案,本发明达到了增强了界面反演的精度和分辨率,并提高了计算效率目的,具有很高的实用价值和推广价值。(The invention discloses a kind of terrestrial heat flow estimation methods based on boundary surface of gravitation and magnetism inverting, this method includes that (S1) carries out inverting to gravity and magnetic data using improved Parker-Oldenburg Interface Inversion algorithm and obtain Curie Point Isothermal Surface and Mohorovich discontinuity, and the inversion algorithm of the Mohorovich discontinuity and Curie Point Isothermal Surface is respectively as follows: Wherein, h 0 For the mean depth of upper and lower interface, Δ h 2 With Δ H 1 For the difference of upper and lower interface and mean depth, s is iteration step length, ρ 0 For the residual density of surface geology medium, ω is circular wavenumber, and F [] is the Fourier transformation of gravity anomaly, and Δ g is gravity anomaly, and G is universal gravitational constant, and a is density with change in depth index；Δ z is magnetic anomaly, μ 0 For magnetic conductivity in vacuum, b is magnetic susceptibility with change in depth index；(S2) underground heat target area is evaluated using Curie Point Isothermal Surface and Mohorovich discontinuity.Through the above scheme, invention achieves enhancing the precision and resolution ratio of Interface Inversion, and computational efficiency purpose is improved, there is very high practical value and promotional value.)

1. a kind of terrestrial heat flow estimation method based on boundary surface of gravitation and magnetism inverting, which comprises the steps of:

(S1) inverting is carried out to gravity and magnetic data using improved Parker-Oldenburg Interface Inversion algorithm and obtains Curie Point Isothermal Surface And Mohorovich discontinuity, the inversion algorithm of the Mohorovich discontinuity and Curie Point Isothermal Surface are respectively as follows:

Wherein, h₀Indicate the mean depth of upper and lower interface, Δ h₂With Δ h₁Indicate the difference of upper and lower interface and mean depth, s is indicated Iteration step length, ρ₀Indicate that the residual density of surface geology medium, ω indicate circular wavenumber, F [] indicates that the Fourier of gravity anomaly becomes It changes, Δ g indicates gravity anomaly, and G indicates universal gravitational constant, and a indicates density with change in depth index；

Wherein, Δ z indicates magnetic anomaly, μ₀Indicate that magnetic conductivity in vacuum, b indicate magnetic susceptibility with change in depth index；

(S2) underground heat target area is evaluated using Curie Point Isothermal Surface and Mohorovich discontinuity.

2. a kind of terrestrial heat flow estimation method based on boundary surface of gravitation and magnetism inverting according to claim 1, which is characterized in that institute Stating the gravity and magnetic data in step (S1) is really quality by changing pole and continuation removes after irregular abnormal data from obtaining matter Measured gravity and magnetic data.

3. a kind of terrestrial heat flow estimation method based on boundary surface of gravitation and magnetism inverting according to claim 2, which is characterized in that institute It states in step (S1) and is included the following steps: using the specific method of improved Parker-Oldenburg Interface Inversion algorithm inverting

(S11) using existing geologic information, geophysical information, drilling data research area lower crust substrate form；

(S12) the irregular magnetic substance magnetic anomaly in upward continuation removal research area's substrate upper caldding layer is utilized；

(S13) binding area base depth and related magnetic susceptibility or density parameter are built with studying area's substrate and the double interfaces of earth's surface Vertical double INTERFACE MODEL bodies；

(S14) the forward modeling gravity-magnetic anomaly of research area's substrate and superstratum is obtained, extraction obtains reflection Curie Point Isothermal Surface or Mohorovich discontinuity rises The deep regional residue gravity-magnetic anomaly of volt；

(S15) using the deep regional residual anomaly obtained in step (S14), using potential field polynomial iterative carry out Curie Point Isothermal Surface or The Inversion Calculation of Mohorovich discontinuity.

4. a kind of terrestrial heat flow estimation method based on boundary surface of gravitation and magnetism inverting according to claim 3, which is characterized in that institute It states in step (S2) by there is transformational relation mathematically using depth of Curie interface and geothermal gradient, together in geothermal exploration When based on Curie Point Isothermal Surface, supplemented by Mohorovich discontinuity, in conjunction with forefathers research data underground heat target area is evaluated.

Technical field

The invention belongs to underground heat neck, geophysics gravimetric prospecting and geophysics magnetic prospecting technical fields, specifically It says, is to be related to a kind of terrestrial heat flow estimation method based on boundary surface of gravitation and magnetism inverting.

Background technique

In recent years, global environmental pollution and energy shortage problem are increasingly severe, are replaced using clean reproducible energy Traditional fossil energy such as traditional petroleum, coal is demand urgently to be resolved.Geothermal energy be exactly it is a kind of cleaning, it is renewable, The widely distributed energy is the ideal target energy.

Geothermal energy resources have a certain depth of burial, thus to geothermal energy resources carry out exploration assessment it is relatively difficult, currently usedly Thermal resource exploitation method is broadly divided into three classes, i.e. geophysical method, geochemical method and boring method.Drilling well is can essence The direct measuring method of the Characteristics of The Distribution of Temperature of underground is really obtained, but boring method, there are well location depth is shallower, valuableness etc. is asked Topic, is not suitable for extensive Reconnaissance Survey of Geothermal Resources Conducted, is only applicable to test process in the target area in geothermal exploration later period.Geophysics class method It is a kind of by quantitative physical method, the method for obtaining the geophysical parameters feature and Geologic Structure Feature of subsurface formations. Herein in the evaluation to geothermal energy resources, the geophysical method used is gravimetric prospecting and magnetic prospecting.

Gravimetric prospecting and magnetic prospecting obtain subsurface rock respectively by the observation to earth internal gravity field and magnetic field Density information and magnetic information.The density and magnetism of rock have sensibility to temperature, therefore can pass through rock density and magnetic Thermo parameters method is heuristically descended in the variation of property.

Corresponding with the heat treatment of underground medium, there is interface Curie Point Isothermal Surfaces for underground.Since pressure acts on, stratum Temperature is increased with the increase of subterranean depth, and on the whole, underground magnetisable material can be increased with temperature demagnetization phenomenon occurs, Cause to will appear the phenomenon that stratum magnetism completely disappears in a certain depth in underground, the temperature front of this depth is referred to as Curie Point Isothermal Surface.The Curie Point Isothermal Surface physical property interface one of important as underground, distribution characteristics directly reflect point of underground temperature field Cloth rule, to tools such as the assessment of geothermal field, the prediction of petroleum resources, the prevention and treatment of earthquake volcanic hazard and primary hydrothermal solution mineral explorations There is great importance.

Corresponding with the variable density of underground medium, there is interface Mohorovich discontinuities for underground.Mohorovich discontinuity is the earth's crust and earth mantle Interface.Mohorovich discontinuity is another important physical property interface of underground, related with the deep factors of ground thermosetting.Study table Bright, Mohorovich discontinuity lifting causes the earth's crust thinning, Upper mantle upwarp, gushes on hot substance, provides source of deep heat for the formation of underground heat.

The fluctuating of Curie Point Isothermal Surface directly reflects the distribution characteristics of underground temperature field, and Mohorovich discontinuity has then affected indirectly ground The distribution of thermal resource, the present invention will with the inverting of Curie Point Isothermal Surface study based on, Mohorovich discontinuity inverting research supplemented by, to geothermal energy resources into Row assessment.

Curie Point Isothermal Surface and the interface undulation of Mohorovich discontinuity can generate remanent magnetization and density, cause the change of earth's surface potential field, symbol The Method And Principle of Parker-Oldenburg Interface Inversion algorithm is closed, therefore Parker-Oldenburg method inverting can be used Curie Point Isothermal Surface and Mohorovich discontinuity.But traditional Parker-Oldenburg method is to carry out Interface Inversion based on single INTERFACE MODEL to grind Study carefully, and what practical geologic body was made of double interfaces or multiple solutions, for this purpose, Wang Wanyin, phase roc et al. are made that corresponding improvement, Inversion result is constrained by using double INTERFACE MODELs.Furthermore physical property (magnetic susceptibility is assumed in traditional Parker-Oldenburg method Or density) parameter is constant, but the variation that the physical parameter in actual formation is complicated in horizontal and vertical upper presence, for this purpose, Researcher carried out improvement, and using varied property model, physical property is expressed as many kinds of function relationship, Feng with change in depth relationship Sharp, the domestic scholars such as Chai Yupu, Jia Jijun will introduce improved Parker- with the varied property parameter of change in depth In Oldenburg method.Improved Parker-Oldenburg inversion algorithm iteration convergence is unstable, although can add Filter makes its convergence, but the use of filter can damage data high-frequency information, reduce the precision of inversion result.Therefore such as What more accurately determines the problem of inverting underground heat target area is those skilled in the art's urgent need to resolve.

Summary of the invention

In order to overcome above-mentioned deficiency in the prior art, the present invention provides a kind of big underground heat based on boundary surface of gravitation and magnetism inverting Flow estimation method, can more accurately inverting Mohorovich discontinuity and Curie Point Isothermal Surface, more perfect evaluation is made to underground heat target area.

To achieve the goals above, The technical solution adopted by the invention is as follows:

A kind of terrestrial heat flow estimation method based on boundary surface of gravitation and magnetism inverting, includes the following steps:

(S1) carrying out inverting to gravity and magnetic data using improved Parker-Oldenburg Interface Inversion algorithm must settle in a foreign country or in another city The inversion algorithm of the inside and Mohorovich discontinuity, the Mohorovich discontinuity and Curie Point Isothermal Surface is respectively as follows:

Wherein, h₀Indicate the mean depth of upper and lower interface, Δ h₂With Δ h₁Indicate the difference of upper and lower interface and mean depth, S indicates iteration step length, ρ₀Indicate that the residual density of surface geology medium, ω indicate circular wavenumber, F [] indicates Fu of gravity anomaly In leaf transformation, Δ g indicate gravity anomaly, G indicate universal gravitational constant, a indicate density with change in depth index；

Wherein, Δ z indicates magnetic anomaly, μ₀Indicate that magnetic conductivity in vacuum, b indicate magnetic susceptibility with change in depth index；

(S2) underground heat target area is evaluated using Curie Point Isothermal Surface and Mohorovich discontinuity.

Further, the gravity and magnetic data in the step (S1) is really that quality is irregular by changing pole and continuation removal From obtaining high-quality gravity and magnetic data after abnormal data.

Further, improved Parker-Oldenburg Interface Inversion algorithm inverting is utilized in the step (S1) Specific method includes the following steps:

(S11) using existing geologic information, geophysical information, drilling data research area lower crust substrate shape State；

(S12) the irregular magnetic substance magnetic anomaly in upward continuation removal research area's substrate upper caldding layer is utilized；

(S13) binding area base depth and related magnetic susceptibility or density parameter, to study area's substrate and the double boundaries of earth's surface Double INTERFACE MODEL bodies are established in face；

(S14) the forward modeling gravity-magnetic anomaly of research area's substrate and superstratum is obtained, extraction obtains reflection Curie Point Isothermal Surface or Moho The deep regional residue gravity-magnetic anomaly that face rises and falls；

(S15) it using the deep regional residual anomaly obtained in step (S14), is occupied using potential field polynomial iterative The Inversion Calculation of the inside or Mohorovich discontinuity.

Specifically, by the way that in geothermal exploration, using depth of Curie interface, there are numbers with geothermal gradient in the step (S4) Transformational relation on:

Wherein, D indicates that depth of Curie interface, Δ T indicate the difference of Curie-point temperature (known) and surface temperature,It indicates Average geothermal gradient can be according to Curie Point Isothermal Surface according to formula (2) it is found that depth of Curie interface and geothermal gradient have inversely prroportional relationship Depth learns the variation of geothermal gradient, understands geothermal distribution, while based on Curie Point Isothermal Surface, supplemented by Mohorovich discontinuity, grinding in conjunction with forefathers The data studied carefully evaluates underground heat target area.

The present invention is by improving after obtaining improvement traditional Parker-Oldenburg Interface Inversion algorithm Parker-Oldenburg Interface Inversion algorithm, specific derivation process are as follows:

Forward Formula is carried out based on the gravity anomaly of dual density interface and vertical variable density in gravimetric prospecting:

It is firstly introduced into variable density parameter, the relationship of density and depth can be assumed as the function of an exponential form:

ρ=ρ₀e^aζ (4)

Wherein, ρ indicates residual density, ρ₀Indicate that the residual density of surface geology medium, a indicate density with change in depth Index, ζ indicate depth, e^aζIndicate density with the change rate of depth；

The bouguer gravity anomaly Δ that a certain density body that underground is at the position (ξ, η, ζ) generates at earth's surface (x, y, 0) G (x, y, 0) are as follows:

Wherein, G indicates universal gravitational constant, and ζ indicates depth (i.e. longitudinal coordinate), and ξ, η then indicate corresponding x, y-coordinate This is reflected in the projected position on ground；

Fourier transformation is carried out to formula (5) to obtain:

Wherein, u and v respectively indicates the wave number on the direction x and y, and i indicates imaginary unit, and D indicates depth of Curie interface, ∫ ∫_D It indicates the limit of integration in triple integral, represents this horizontal cross face, formula (6) convert:

If interface is h in geologic body₂(x, y), lower interface are h₁(x, y), then the mean depth of upper and lower interface is h₀, on The difference of lower interface and mean depth is respectively Δ h₂With Δ h₁, enableThen formula (7) is represented by as follows Form:

By e^(a-ω)ζThe Taylor expansion at ζ=0, and ζ is integrated:

Formula (9) is recombinated, following formula is obtained:

Formula (9) the as Parker-Oldenburg gravity anomaly based on dual density interface and vertical variable density is just Drill formula.

Inversion formula is carried out to the lower interface under dual density INTERFACE MODEL variable density constraint in gravimetric prospecting:

Interface Inversion formula is obtained for iterative calculation, it is first assumed that the gravity anomaly initial value of a certain depth in underground is set For Δ g (x, y, z)⁽¹⁾=Δ g (x, y, 0) show that the exceptional value of earth's surface calculates initial value using formula (10) are as follows:

In conjunction with potential field iterative calculation method, obtain:

Δ g (x, y, z)⁽²⁾=Δ g (x, y, z)⁽¹⁾+ s (Δ g (x, y, 0)-Δ g (x, y, 0)⁽¹⁾) (12)

Wherein, s indicates iteration step length, and so on, the result after n times iterate to calculate:

Δ g (x, y, z)⁽ⁿ⁺¹⁾=Δ g (x, y, z)⁽ⁿ⁾+ s (Δ g (x, y, 0)-Δ g (x, y, 0)⁽ⁿ⁾) (13)

When | Δ g (x, y, 0)-Δ g (x, y, 0)⁽ⁿ⁾When |≤ε, it is believed that Δ g (x, y, z)⁽ⁿ⁺¹⁾≈Δg(x,y,z)⁽ⁿ⁾, to public affairs Formula (13) derives:

Depth results are extracted to obtain:

Formula (15) is the lower interface inversion formula under the constraint of dual density INTERFACE MODEL variable density.

Magnetic anomaly Forward Formula based on double magnetic interfaces and vertical metamagnetism:

It is firstly introduced into variable magnetization rate parameter, the relationship of magnetic susceptibility and depth can be assumed as the function of an exponential form:

M=M₀e^aζ (16)

Wherein, M₀Indicate the remanent magnetization of surface geology medium, a indicates magnetic susceptibility with change in depth index, ζ table Show depth, e^aζIndicate magnetic susceptibility with the change rate of depth；

It can be obtained by frequency domain Poisson formula:

Wherein, F [] indicates the Fourier transformation of magnetic anomaly, U_mIndicate magnetic potential, V indicates that gravity position, M indicate remanent magnetization Intensity, ω indicate circular wavenumber；G indicates universal gravitational constant, and ρ indicates residual density；

Formula 16 is substituted into 17, following formula is obtained:

The wave spectrum of magnetic anomaly can indicate are as follows:

F [Δ z]=- μ₀ωF[U_m] (19)

Wherein, Δ z indicates magnetic anomaly, μ₀Indicate magnetic conductivity in vacuum；

The wave spectrum of gravity anomaly can be expressed as form:

F [Δ g]=ω F [V] (20)

Aggregative formula (17) (18) (19) (20) can obtain:

Formula (8) substitution formula (21) can be obtained:

Wherein, i indicates imaginary unit, by e^(a-ω)ζIn ζ=h₀Locate Taylor expansion, and ζ integrated:

Formula (23) is recombinated, following formula is obtained:

Formula (23) is perpendicular magnetization rate and double magnetic interface Model Condition lower frequency domain forward modeling magnetic anomaly formula.

Lower interface inversion formula under double magnetic interface model variable magnetization rate constraints:

Interface Inversion formula is obtained for iterative calculation, it is first assumed that the gravity anomaly initial value of a certain depth in underground is set For Δ z (x, y, z)⁽¹⁾=Δ z (x, y, 0) show that the exceptional value of earth's surface calculates initial value using formula (23) are as follows:

In conjunction with potential field iterative calculation method, obtain:

Δz(x,y,z)⁽²⁾=Δ z (x, y, z)⁽¹⁾+s(Δz(x,y,0)-Δz(x,y,0)⁽¹⁾) (26)

Wherein, s indicates iteration step length, and so on, the result after n times iterate to calculate:

Δ z (x, y, z)⁽ⁿ⁺¹⁾=Δ z (x, y, z)⁽ⁿ⁾+ s (Δ z (x, y, 0)-Δ z (x, y, 0)⁽ⁿ⁾) (27)

When | Δ z (x, y, 0)-Δ z (x, y, 0)⁽ⁿ⁾When |≤ε, it is believed that Δ z (x, y, z)⁽ⁿ⁺¹⁾≈Δz(x,y,z)⁽ⁿ⁾, to public affairs Formula (26) derivation can obtain:

Depth results are extracted to obtain:

Formula (29) (i.e. formula (1)) is the lower interface inversion formula under double magnetic interface model variable magnetization rate constraints.

Compared with prior art, the invention has the following advantages:

(1) present invention improves traditional Parker-Oldenburg method, and traditional single INTERFACE MODEL is improved For double INTERFACE MODELs (i.e. comprising Curie Point Isothermal Surface and Mohorovich discontinuity), double Interface Inversions are realized, it is opposite with tradition by double interface detections Single interface detection, test data is more accurate, and improved Parker-Oldenburg method is when testing, draws The varied property restriction on the parameters for having entered more to meet underground medium physical property obtains to obtain the data for meeting practical boundary surface of gravitation and magnetism The case where therefore data are more in line with underground medium, and calculated using potential field polynomial iterative, the data result obtained is more Accurately, to improve the precision and reliability of inversion result, theories integration is provided for terrestrial heat flow estimation.

(2) new iterative algorithm is introduced into Parker-Oldenburg Interface Inversion formula by the present invention, improves inverting As a result precision, while new iterative algorithm is guaranteeing the convergent premise of algorithm iteration process, avoids making for bandpass filter With remaining the effective high-frequency information of initial data in refutation process, improve the accuracy and integrality of inversion result. And the algorithm computational efficiency proposed in the present invention improves four times or more relative to traditional technical approach on calculation rate.

Detailed description of the invention

Fig. 1 is system structure flow chart of the invention.

Fig. 2 is result figure after magnetic anomaly normalizing pole continuation in republicanism basin of the present invention.

Fig. 3 is republicanism basin of the present invention base depth map.

Fig. 4 is republicanism basin of the present invention substrate and superstratum forward modeling magnetic anomaly map.

Fig. 5 is the region residue magnetic anomaly map that republicanism basin of the present invention reflects that Curie Point Isothermal Surface rises and falls.

Fig. 6 is Curie Point Isothermal Surface inversion result figure in republicanism basin of the present invention.

Fig. 7 is result figure after republicanism basin gravity anomaly upward continuation of the present invention.

Fig. 8 is republicanism basin of the present invention substrate and superstratum forward modeling gravity anomaly map.

Fig. 9 is the region residual gravity anomaly figure that republicanism basin of the present invention reflects that Curie Point Isothermal Surface rises and falls.

Figure 10 is Mohorovich discontinuity inversion result figure in republicanism basin of the present invention.

Specific embodiment

Present invention will be further explained below with reference to the attached drawings and examples, and embodiments of the present invention include but unlimited In the following example.