阅读说明：本技术 一种基于重磁界面反演的大地热流估计方法 (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.)

一种基于重磁界面反演的大地热流估计方法

技术领域

本发明属于地热领、地球物理重力勘探和地球物理磁法勘探技术领域，具 体地讲，是涉及一种基于重磁界面反演的大地热流估计方法。

背景技术

近年来，全球环境污染和能源短缺问题越来越严重，使用清洁可再生能源 来代替传统的石油、煤炭等传统的化石能源，是亟待解决的需求。地热能就是 一种清洁、可再生、分布广泛的能源，是理想的目标能源。

地热资源具有一定埋深，因此对地热资源进行勘探评估比较困难，目前常 用的地热资源勘探方法主要分为三类，即地球物理方法，地球化学方法和钻井 方法。钻井是能精确获取地下的温度分布特征的直接测量方法，但钻井方法存 在井位深度较浅，昂贵等问题，不适合大规模地热普查，仅适用于于地热勘探 后期的靶区试验过程。地球物理类方法是一种通过定量的物理方法，获取地下 地层的地球物理参数特征和地质构造特征的方法。本文在对地热资源的评价中， 采用的地球物理方法为重力勘探和磁法勘探。

重力勘探和磁法勘探分别通过对地球内部重力场和磁场的观测，获取地下 岩石的密度信息和磁性信息。岩石的密度和磁性对温度具有敏感性，因此能通 过岩石密度和磁性的变化推断地下温度场分布。

与地下介质的磁性变化相对应的，地下存在着界面居里面。由于压力作用， 地层温度随着地下深度的增加而升高，整体上，地下磁性物质会随着温度升高 出现退磁现象，导致在地下某一深度处会出现地层磁性完全消失的现象，这一 深度处的温度界面被称为居里面。居里面作为地下重要的物性分界面之一，其 分布特征直接反映了地下温度场的分布规律，对地热田的评估、油气资源的预 测、地震火山灾害防治及原生热液矿产勘查等具有重要指导意义。

与地下介质的密度变化相对应的，地下存在着界面莫霍面。莫霍面是地壳 和地幔的分界面。莫霍面是地下另一个重要的物性分界面，与地热形成的深部 因素有关。研究表明，莫霍面抬升导致地壳变薄，上地幔***，热物质上涌， 为地热的形成提供了深部热源。

居里面的起伏直接反映了地下温度场的分布特征，莫霍面则间接地影响了 了地热资源的分布，本发明将以居里面的反演研究为主，莫霍面反演研究为辅， 对地热资源进行评估。

居里面和莫霍面的界面起伏均会产生剩余磁性和密度，引起地表位场的改 变，符合Parker-Oldenburg界面反演算法的方法原理，因此可以使用 Parker-Oldenburg法反演居里面和莫霍面。但是，传统的Parker-Oldenburg法是 基于单界面模型开展界面反演研究，而实际地质体是由双界面或多界面组成的， 为此，王万银，相鹏等人做出了相应改进，通过采用双界面模型约束反演结果。 此外传统的Parker-Oldenburg法中假设物性(磁化率或密度)参数为常数，但是实 际地层中的物性参数在横向和纵向上存在复杂的变化，为此，研究人员进行过 改进，采用变物性模型，将物性随深度变化关系表示成多种函数关系，冯锐， 柴玉璞，贾继军等国内学者将随深度变化的变物性参数引入到改进后的 Parker-Oldenburg方法中。改进后的Parker-Oldenburg反演算法迭代收敛不稳定， 虽然可以添加滤波器来使其收敛，但滤波器的使用会损害数据高频信息，降低 反演结果的精度。因此如何更加准确地确定反演地热靶区是本领域技术人员亟 需解决的问题。

发明内容

为了克服现有技术中的上述不足，本发明提供一种基于重磁界面反演的大 地热流估计方法，能够更准确地反演莫霍面和居里面，对地热靶区做出更完善 的评价。

为了实现上述目的，本发明采用的技术方案如下：

一种基于重磁界面反演的大地热流估计方法，包括如下步骤：

(S1)利用改进后的Parker-Oldenburg界面反演算法对重磁数据进行反演得 出居里面和莫霍面，所述莫霍面和居里面的反演算法分别为：

其中，h₀表示上下界面的平均深度，Δh₂和Δh₁表示上下界面与平均深度的 差值，s表示迭代步长，ρ₀表示地表地质介质的剩余密度，ω表示圆波数，F[] 表示重力异常的傅里叶变换，Δg表示重力异常，G表示万有引力常数，a表示密 度随深度变化指数；

其中，Δz表示磁异常，μ₀表示真空中磁导率，b表示磁化率随深度变化指 数；

(S2)利用居里面和莫霍面对地热靶区进行评价。

进一步，所述步骤(S1)中的重磁数据实际是质量通过化极和延拓去除不 规则的异常数据后从得到质量好的重磁数据。

进一步地，所述步骤(S1)中利用改进的Parker-Oldenburg界面反演算法反 演的具体方法包括如下步骤：

(S11)利用已有的地质资料、地球物理资料、钻井资料研究区地下地壳基 底的形态；

(S12)利用向上延拓去除研究区基底上覆盖层中的不规则磁性体磁异常；

(S13)结合研究区基底深度和相关磁化率或密度参数，以研究区基底和地 表双界面建立双界面模型体；

(S14)获取研究区基底及上覆地层的正演重磁异常，提取得到反映居里面 或莫霍面起伏的深部区域剩余重磁异常；

(S15)利用步骤(S14)中获取的深部区域剩余异常，采用位场多项式迭代 进行居里面或莫霍面的反演计算。

具体地，所述步骤(S4)中通过在地热勘探中，利用居里面深度与地温梯 度存在数学上的转换关系：

其中，D表示居里面深度，ΔT表示居里点温度(已知)与地表温度的差值， 表示平均地温梯度，根据公式(2)可知，居里面深度和地温梯度具有反比例 关系，可根据居里面深度得知地温梯度的变化，了解地热分布，同时以居里面 为主，莫霍面为辅，再结合前人研究的资料对地热靶区进行评价。

本发明通过对传统的Parker-Oldenburg界面反演算法进行改进得到改进后Parker-Oldenburg界面反演算法，其具体推导过程如下：

在重力勘探中基于双密度界面和垂向变密度的重力异常进行正演公式：

首先引入变密度参数，密度和深度的关系可以假定为一个指数形式的函数：

ρ＝ρ₀e^aζ (4)

其中，ρ表示剩余密度，ρ₀表示地表地质介质的剩余密度，a表示密度随深 度变化指数，ζ表示深度，e^aζ表示密度随深度的变化率；

地下处于(ξ,η,ζ)位置处的某一密度体在地表(x,y,0)处产生的布格重力异常 Δg(x,y,0)为：

其中，G表示万有引力常数，ζ表示深度(即纵向坐标)，ξ、η则表示对应 的x，y坐标反映了该点在地面的投影位置；

对公式(5)进行傅里叶变换得到：

其中，u和v分别表示x和y方向上的波数，i表示虚数单位，D表示居里 面深度，∫∫_D表示三重积分中的积分限，代表这水平横向面，对公式(6)进行变换 得：

如果地质体上界面为h₂(x,y)，下界面为h₁(x,y)，则上下界面的平均深度为 h₀，上下界面与平均深度的差值分别为Δh₂和Δh₁，令则公式(7) 可表示为如下形式：

将e^(a-ω)ζ在ζ＝0处泰勒展开，并对ζ积分得:

对公式(9)进行重组，得到如下公式：

公式(9)即为基于双密度界面和垂向密度变化的Parker-Oldenburg重力异常 正演公式。

对重力勘探中双密度界面模型变密度约束下的下界面进行反演公式：

为迭代计算得出界面反演公式，首先假定地下某一深度处的重力异常初始 值设为Δg(x,y,z)⁽¹⁾＝Δg(x,y,0)，利用公式(10)得出地表的异常值计算初值为：

结合位场迭代计算方法，得：

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

其中，s表示迭代步长，以此类推，经过n次迭代计算后的结果：

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

当|Δg(x,y,0)-Δg(x,y,0)⁽ⁿ⁾|≤ε时，认为Δg(x,y,z)⁽ⁿ⁺¹⁾≈Δg(x,y,z)⁽ⁿ⁾，对公式(13)推导得：

提取深度结果得到：

公式(15)即为双密度界面模型变密度约束下的下界面反演公式。

基于双磁性界面和垂向变磁性的磁异常正演公式：

首先引入变磁化率参数，磁化率和深度的关系可以假定为一个指数形式的 函数：

M＝M₀e^aζ (16)

其中，M₀表示地表地质介质的剩余磁化强度，a表示磁化率随深度变化指 数，ζ表示深度，e^aζ表示磁化率随深度的变化率；

由频率域泊松公式可得:

其中，F[]表示磁异常的傅里叶变换，U_m表示磁位，V表示重力位，M表示 剩余磁化强度，ω表示圆波数；G表示万有引力常数，ρ表示剩余密度；

将公式16代入17中，获得下式：

磁异常的波谱可以表示为：

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

其中，Δz表示磁异常，μ₀表示真空中磁导率；

重力异常的波谱可以表示为如下形式：

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

综合公式(17)(18)(19)(20)可得：

将公式(8)代入公式(21)即可得:

其中，i表示虚数单位，将e^(a-ω)ζ在ζ＝h₀处泰勒展开，并对ζ积分得:

对公式(23)进行重组，得到如下公式：

公式(23)即为垂向磁化率和双磁性界面模型条件下频率域正演磁异常公式。

双磁性界面模型变磁化率约束下的下界面反演公式：

为迭代计算得出界面反演公式，首先假定地下某一深度处的重力异常初始 值设为Δz(x,y,z)⁽¹⁾＝Δz(x,y,0)，利用公式(23)得出地表的异常值计算初值为：

结合位场迭代计算方法，得：

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

其中，s表示迭代步长，以此类推，经过n次迭代计算后的结果：

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

当|Δz(x,y,0)-Δz(x,y,0)⁽ⁿ⁾|≤ε时，认为Δz(x,y,z)⁽ⁿ⁺¹⁾≈Δz(x,y,z)⁽ⁿ⁾，对公式 (26)推导可得：

提取深度结果得：

公式(29)(即公式(1))是双磁性界面模型变磁化率约束下的下界面反演 公式。

与现有技术相比，本发明具有以下有益效果：

(1)本发明对传统的Parker-Oldenburg方法进行改进，将传统的单界面 模型改进为双界面模型(即包含居里面和莫霍面)，实现了双界面反演，通过 双界面测试相对与传统的单界面测试，测试数据更精准，并且改进后的 Parker-Oldenburg方法在进行测试时，引入了更符合地下介质物性的变物性参 数约束，从而得出符合实际重磁界面的数据，得出的数据因此更加符合地下介 质的情况，并且采用位场多项式迭代计算，得出的数据结果更加精确，从而提 高了反演结果的精度和可靠性，为大地热流估计提供了理论支持。

(2)本发明将新迭代算法引入到Parker-Oldenburg界面反演公式中，提高 了反演结果的精度，同时新迭代算法在保证算法迭代过程收敛的前提，避免了 带通滤波器的使用，在反演过程中保留了原始数据有效的高频信息，提高了反 演结果的准确性和完整性。且本发明中提出的算法计算效率相对于传统的技术 方式在计算率上提高四倍及以上。

附图说明

图1为本发明的系统结构流程图。

图2为本发明共和盆地磁异常化极延拓后结果图。

图3为本发明共和盆地基地深度图。

图4为本发明共和盆地基底及上覆地层正演磁异常图。

图5为本发明共和盆地反映居里面起伏的区域剩余磁异常图。

图6为本发明共和盆地居里面反演结果图。

图7为本发明共和盆地重力异常向上延拓后结果图。

图8为本发明共和盆地基底及上覆地层正演重力异常图。

图9为本发明共和盆地反映居里面起伏的区域剩余重力异常图。

图10为本发明共和盆地莫霍面反演结果图。

具体实施方式

下面结合附图和实施例对本发明作进一步说明，本发明的实施方式包括但 不限于下列实施例。

实施例

如图1所示，一种基于重磁界面反演的大地热流估计方法，包括如下步骤：

(S1)利用改进后的Parker-Oldenburg界面反演算法对重磁数据进行反演得 出居里面和莫霍面，所述莫霍面和居里面的反演算法分别为：

其中，h₀表示上下界面的平均深度，Δh₂和Δh₁表示上下界面与平均深度的 差值，s表示迭代步长，ρ₀表示地表地质介质的剩余密度，ω表示圆波数，F[] 表示重力异常的傅里叶变换，Δg表示重力异常，G表示万有引力常数，a表示密 度随深度变化指数；

其中，Δz表示磁异常，μ₀表示真空中磁导率，b表示磁化率随深度变化指 数；

（S2）利用居里面和莫霍面对地热靶区进行评价。

进一步，所述步骤（S1）中的重磁数据实际是质量通过化极和延拓去除不 规则的异常数据后从得到质量好的重磁数据。

进一步地，所述步骤（S1）中利用改进的Parker-Oldenburg界面反演算法反 演的具体方法包括如下步骤：

（S11）利用已有的地质资料、地球物理资料、钻井资料研究区地下地壳基 底的形态；

（S12）利用向上延拓去除研究区基底上覆盖层中的不规则磁性体磁异常；

（S13）结合研究区基底深度和相关磁化率或密度参数，以研究区基底和地 表双界面建立双界面模型体；

（S14）获取研究区基底及上覆地层的正演重磁异常，提取得到反映居里面 或莫霍面起伏的深部区域剩余重磁异常；

（S15）利用步骤（S14）中获取的深部区域剩余异常，采用位场多项式迭代 进行居里面或莫霍面的反演计算。

具体地，所述步骤（S4）中通过在地热勘探中，利用居里面深度与地温梯 度存在数学上的转换关系：

其中，D表示居里面深度，ΔT表示居里点温度（已知）与地表温度的差值， 表示平均地温梯度，根据公式（2）可知，居里面深度和地温梯度具有反比例 关系，可根据居里面深度得知地温梯度的变化，了解地热分布，同时以居里面 为主，莫霍面为辅，再结合前人研究的资料对地热靶区进行评价。

为了测试本发明实际应用效果，以共和盆地地热资源勘探为例，利用本发 明中提出的反演算法，对共和盆地的居里面和莫霍面进行反演，并进行地热评 价

先进行居里面的反演，共和盆地航磁异常值在-120—180nT之间波动，磁异 常主要分布方向为NE向和NW向。图中盆地及周边地区磁异常展布方向为NE 向，盆地内部磁异常展布为NW向，呈NW向高ES向低，由NW向至SE向呈 逐渐降低的态势。图2是对初始磁异常进行化极和延拓后的结果，图3是根据 地震数据和钻井资料计算出的基地深度图。根据基地深度，利用公式(22)计 算出基底及上覆地层的正演磁异常，正演磁异常如图4所示。再从图2的总磁 异常中减去图4中基底及上覆地层的磁异常，得到如图5所示的反映居里面起 伏的剩余磁异常。利用公式(1)即可反演得出居里面深度，其中变磁化率参数 通过井深温度分布资料获得，以青海省地层的平均磁化率1.3847×10-3SI为起 算点，得到基底与居里面的之间的磁化率随深度的变化关系为κ＝1.97×10^-5e^-0.1z， 反演结果如图6所示。

接着进行莫霍面的反演，共和盆地物性资料显示，由钻井资料获得了共和 盆地剩余岩石密度参数为0.6g/cm³。图7是对初始重力异常进行了向上延拓处理， 去除了浅部的一些不规则重力异常后的结果。根据图3所示基底深度，利用公 式(9)计算出基底及上覆地层的正演重力异常，正演重力异常如图8所示。再 从图7的重力异常中减去图8中基底及上覆地层的重力异常，得到如图9所示 的反映莫霍面起伏的剩余重力异常。利用公式(15)即可反演得出莫霍面深度， 反演结果如图10所示。

最后根据获得的居里面和莫霍面的反演结果，对地热靶区做出评价。根据 公式(3)，地温梯度和居里面具有反比关系。通过对共和盆地居里面和莫霍面 反演结果进行分析，可以推断盆地及周边区域居里面和莫霍面形变主要受NW 向展布的断裂构造影响；盆地内居里面和莫霍面形变主要受NE向展布的断裂构 造影响控制。图中显示共和盆地中部居里面和莫霍面反演结果均较浅，推测莫 霍面和居里面***是共和盆地地温梯度和大地热流值高的深层次原因。居里面 和莫霍面上隆现象为共和盆地地热的形成提供了有利的深部热源条件；此外， 盆地内发育两条深大断裂，可以成为地下热量向上传递的良好通道。盆地内具 有较厚的沉积盖层，可以很好的保存热量，由此推断共和盆地具有良好的地热 开发远景。

通过上述方法，本发明能很好地提高重磁数据反演中居里面和莫霍面的反 演精度，提高计算效率，以此更准确地圈定地热靶区。因此，本发明具有突出 的实质性的特点和进步。

上述实施例仅为本发明的优选实施例，并非对本发明保护范围的限制，但 凡采用本发明的设计原理，以及在此基础上进行非创造性劳动而做出的变化， 均应属于本发明的保护范围之内。