专利摘要

本发明公开了一种基于重磁界面反演的大地热流估计方法，该方法包括(S1)利用改进后的Parker‑Oldenburg界面反演算法对重磁数据进行反演得出居里面和莫霍面，所述莫霍面和居里面的反演算法分别为：其中，h0为上下界面的平均深度，Δh2和Δh1为上下界面与平均深度的差值，s为迭代步长，ρ0为地表地质介质的剩余密度，ω为圆波数，F[]为重力异常的傅里叶变换，Δg为重力异常，G为万有引力常数，a为密度随深度变化指数；Δz为磁异常，μ0为真空中磁导率，b为磁化率随深度变化指数；(S2)利用居里面和莫霍面对地热靶区进行评价。通过上述方案，本发明达到了增强了界面反演的精度和分辨率，并提高了计算效率目的，具有很高的实用价值和推广价值。

权利要求

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

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

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

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

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

2.根据权利要求1所述的一种基于重磁界面反演的大地热流估计方法，其特征在于，所述步骤(S1)中的重磁数据实际是通过化极和延拓去除不规则的异常数据后从而得到质量好的重磁数据。

3.根据权利要求2所述的一种基于重磁界面反演的大地热流估计方法，其特征在于，所述步骤(S1)中利用改进的Parker-Oldenburg界面反演算法反演的具体方法包括如下步骤：

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

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

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

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

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

4.根据权利要求3所述的一种基于重磁界面反演的大地热流估计方法，其特征在于，所述步骤(S2)中通过在地热勘探中，利用居里面深度与地温梯度存在数学上的转换关系，同时以居里面为主，莫霍面为辅，再结合前人研究的资料对地热靶区进行评价。

说明书

技术领域

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

背景技术

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

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

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

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

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

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

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

发明内容

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

ρ＝ρ0eaζ(4)

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

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

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

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

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

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

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

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

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

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

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

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

Δg(x,y,z)(2)＝Δg(x,y,z)(1)+s(Δg(x,y,0)-Δg(x,y,0)(1))(12)

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

Δg(x,y,z)(n+1)＝Δg(x,y,z)(n)+s(Δg(x,y,0)-Δg(x,y,0)(n))(13)

当|Δg(x,y,0)-Δg(x,y,0)(n)|≤ε时，认为Δg(x,y,z)(n+1)≈Δg(x,y,z)(n)，对公式(13)推导得：

提取深度结果得到：

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

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

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

M＝M0ebζ(16)

其中，M0表示地表地质介质的剩余磁化强度，b表示磁化率随深度变化指数，ζ表示深度，ebζ表示磁化率随深度的变化率；

由频率域泊松公式可得:

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

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

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

F[Δz]＝-μ0ωF[Um] (19)

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

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

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

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

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

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

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

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

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

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

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

Δz(x,y,z)(2)＝Δz(x,y,z)(1)+s(Δz(x,y,0)-Δz(x,y,0)(1))(26)

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

Δz(x,y,z)(n+1)＝Δz(x,y,z)(n)+s(Δz(x,y,0)-Δz(x,y,0)(n)) (27)

当|Δz(x,y,0)-Δz(x,y,0)(n)|≤ε时，认为Δz(x,y,z)(n+1)≈Δz(x,y,z)(n)，对公式(26)推导可得：

提取深度结果得：

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

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

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

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

附图说明

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

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

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

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

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

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

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

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

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

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

具体实施方式

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

实施例

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

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

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

其中，Δz表示磁异常，μ0表示真空中磁导率，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/cm3。图7是对初始重力异常进行了向上延拓处理，去除了浅部的一些不规则重力异常后的结果。根据图3所示基底深度，利用公式(9)计算出基底及上覆地层的正演重力异常，正演重力异常如图8所示。再从图7的重力异常中减去图8中基底及上覆地层的重力异常，得到如图9所示的反映莫霍面起伏的剩余重力异常。利用公式(15)即可反演得出莫霍面深度，反演结果如图10所示。

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

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

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

一种基于重磁界面反演的大地热流估计方法专利购买费用说明

Q：办理专利转让的流程及所需资料

A：专利权人变更需要办理著录项目变更手续，有代理机构的，变更手续应当由代理机构办理。

1：专利变更应当使用专利局统一制作的“著录项目变更申报书”提出。

2：按规定缴纳著录项目变更手续费。

3：同时提交相关证明文件原件。

4：专利权转移的，变更后的专利权人委托新专利代理机构的，应当提交变更后的全体专利申请人签字或者盖章的委托书。

Q:专利著录项目变更费用如何缴交

A：（1）直接到国家知识产权局受理大厅收费窗口缴纳,（2）通过代办处缴纳，（3）通过邮局或者银行汇款,更多缴纳方式

Q：专利转让变更，多久能出结果

A：著录项目变更请求书递交后，一般1-2个月左右就会收到通知，国家知识产权局会下达《转让手续合格通知书》。