谷露地热田热结构初探及高温异常区预测

王闻文，巴桑元旦，吴儒杰

1 核工业北京地质研究院，北京 100029

2 西藏自治区土地矿权交易和资源储量评审中心，西藏 拉萨 850011

温度作为7 个基本物理量之一，是极其重要的物理参数，地球内部的温度结构不仅影响大陆变形和构造演化过程［1-2］，而且决定了浅地表地热资源的分布。地热能作为绿色清洁能源，尤其是高温地热，具有稳定连续输出的发电优势，近年来在全球范围内得到了广泛的利用［3-4］。地热田地下1～2 km 深度内，地下介质的温度及其分布规律直接决定了高温热储的界定、利用与开发［5］。因此，揭示地下温度结构特征对于地热田的研究、勘探和开发至关重要。

目前，地球内部温度数据的获取一般包括直接和间接两种方法［6］。直接的方法是钻孔测温法，即通过测温仪来对井液的温度进行测量，当钻孔经过相对较长的静井恢复时间后，井液温度与地层温度达到热平衡状态［7］，此时的地温数据，可以较为真实地代表对应地层的温度。钻孔测温也是地热田的勘查过程中使用最为广泛且最为首要的方法。然而，由于受到钻孔施工的成本、时间周期和技术难度等因素的限制，钻孔测温法往往难于批量地投入生产。

相比于直接法，采用间接方法来估算地温数据往往具有更低的经济成本、更好的空间分布连续性和更大的地层深度。常用的地温间接测量法，主要可分为两类［8］：1）热传导建模［9-11］； 2）间接地温计［12］。第1 类方法构建温度模型，首先需要确定传热机制（即确定传热是通过传导、对流还是两者结合），其次是定义边界条件，最后是准确定义各种介质的热物性参数；其中，热物性参数和边界条件（特别是下边界条件）只能进行粗略估算。第2 类方法是在地表收集地质、地球化学、同位素或气体成分数据的基础上加以分析，如典型的地热温标法［13］，其虽然可以评估某些特征深度的温度，但是却很难获得地下空间区域范围内温度场的分布规律。

近年来的研究认为［6，8］，要认识地温场的空间分布特征，可以通过建立温度和与其相关的岩石物性参数之间的经验公式来实现。岩石的物性参数中，电阻率是对温度反应最敏感的参数，这是因为地热系统通常含有含盐流体，而水热蚀变过程导致系统发育的岩石电阻率发生普遍变化。具体而言，含盐流体在水热蚀变或高温作用下，往往会导致地热系统的电阻率降低，与围岩形成明显的电性差异，易于被电阻率测深法成像［14］。由此，基于电阻率数据约束的温度预测法，被广泛应用于地热资源的勘探［15-19］。

本文充分利用了电阻率法在温度预测中的优势以及音频大地电磁测深法（Audiofrequency Magnetotellurics，简称AMT）在探测地热田地电结构中的有效性，采用Campbell、Keller 和Chermensky 等建立的经验公式［20-22］来定义岩土电阻率与温度的关系，并基于谷露地热田的钻孔测温数据，通过回归分析建立了温度与深度的线性函数关系，随即对温度相关系数根据不同深度进行了分层计算，构建了不同空间坐标下的温度预测模型，通过逐层计算地下温度，建立了谷露地热田L1、L2 两条AMT 反演电阻率剖面的温度结构模型。在此基础上，通过两个验证钻孔YZ01、YZ02 处的温度预测值与钻孔测温实际值的对比，结合谷露地热田地质特征及前人研究成果，初步证实了采用温度估算方法的准确性和开展高温异常区预测的有效性。将电阻率勘探方法及钻孔测温法在地热田勘查过程中深度融合，形成具有优势互补的特点，并充分挖掘两者之间的内在联系，对于深入揭示类似于谷露地热地质背景的其他地热田的热结构具有重要的启示意义。

1 地质背景

西藏谷露地热田位于亚东—谷露裂谷的北部（图1a）、谷露盆地中部，隶属于那曲—尼木地热带，带内分布有续迈、羊易、羊八井、谷露和玉寨等26 个中高温地热显示区［23-24］（图1b），这些地热显示区均出露于山前“串珠”状断陷盆地内，整体受念青唐古拉东南麓活动断裂控制［25］。热田成因受地球深部热动力学过程、地表断裂活动等因素的控制，地热带内在18.3～11.1 Ma 期间发生大规模岩浆侵入，形成了念青唐古拉巨大花岗岩岩基［26］。

图1 那曲—尼木地热带大地构造及地热显示区分布（据参考文献 ［23］ 修改）Fig. 1 Tectonic of Naqu-Nimu geothermal belt and distribution of geothermal field

谷露地热田（图2）两侧为中低山区：西侧为花岗闪长岩体、马里组（J2m）灰岩与砂岩互层的变质岩；东侧地势较为平坦，主要为第四系（Q4）松散堆积物［27-30］。花岗闪长岩体侵入至马里组中，局部出露规模较小的花岗斑岩体和酸性岩脉，侵入岩时代均为新生代［24］，并构成区内第四系沉积地层的基底［31］。

图2 谷露地热田断裂构造与地热形迹（据参考文献 ［24］ 修改）Fig. 2 Structure and geothermal traces of Gulu Geothermal Field

谷露地热田内主要发育EW 向、NE 向和SN向3 组断裂（图2）。EW 向F1-1和F1-3断裂属区域性断裂，其中F1-1断裂全长约5 km，贯穿整个谷露盆地；F1-3断裂从泉华台地北侧通过，同样贯穿盆地，向西一直延伸至念青唐古拉山，在盆地西侧表现为宽缓的断层谷地，南北最宽处约900 m，谷地内大面积分布中上更新统冲洪积地层［24］。NE 向断裂（F2-1～F2-9）出露于热田基底花岗岩中［23］，自西南往东北，连通念青唐古拉山；结合野外地质调查及大地电磁测深反演结果，认为NE 向断裂是谷露地热田地下水补给与热水运移的重要通道［24］。SN 向断裂主要发育3 条（F4-1～F4-3），为区域张性断裂，其中F4-1断裂位于盆地西侧花岗岩中，在地表呈SN向断层陡坎地貌；F4-2断裂在地表特征最为明显，沿桑曲河西岸发育，断层倾向E，倾角介于60°～70°之间，热水活动皆分布于F4-2断裂上盘；F4-3断裂位于桑曲河东岸，长约2 km，走向介于10°～15°之间，在地表断层倾角近陡立，以西倾为主，错断泉胶砾岩层，沿断层热水泉眼线状出露；地表地热显示区集中分布于F4-2和F4-3断裂的夹持区域，认为SN 向断裂F4-2导通深部热源，是深部热流体上移的重要通道［23-24］。

2 研究方法

2.1 电阻率与温度转换关系的建立

电阻率是岩石非常重要的物理性质，已被证明是寻找地热资源时最有用的地球物理参数之一，在地热资源的勘探过程中被广泛利用。如何在地热田建立温度与电阻率之间的转换关系，就显得尤为重要。

温度对电阻率的影响，源于温度会改变岩石本身内在的物理性质，从而对其导电性能产生影响。岩石的电阻率与其物理性质这二者之间的经验公式中，Archie 公式［32］影响深远，在石油测井中被广泛使用。地层因素的Archie公式［33］为：

式（1）中：R0—100%饱含地层水时的地层电阻率，Ω·m；Rw—地层水的电阻率，Ω·m；Φ—孔隙度，无量纲；m—岩石的胶结指数，无量纲；a—与岩石有关的比例系数，无量纲。

实质上，Archie 利用岩石电阻率与孔隙水电阻率的比值来反映地层因素，建立了岩石物性和电性之间的联系［34］。借鉴电阻率比值这一思想，Campbell、Keller 和Cheremensk 等众多学者先后建立了不同的经验公式［20-22］来反映温度变化对岩石电阻率的影响，但总体可以用公式（2） 来表达：

式（2）中：RT— 给定温度下岩石的电阻率，Ω·m；Ts—给定的温度值，℃；αs—Ts温度值下的温度系数，无量纲；Rs—Ts温度值下的岩石标准电阻率，Ω·m。

通过李艳华等［35］对深约1.8 km 的六块砂岩岩心进行的岩石物理实验，可以直观地了解温度对岩石电阻率的影响程度。图3 中展示了6 条呈幂指数衰减的曲线，反映了这六块砂岩岩心的电阻率随着温度升高而呈现降低的变化趋势。此外，Caldwell 等人［36］也早在1986 年研究了东南亚地区地热钻孔岩心的电阻率与温度间的关系，并通过回归分析得到了相关系数的值。Anderson 等人［14］在Caldwell 等人的基础上进一步研究了热储层中岩石电阻率与温度间的关系。通过上述学者的研究表明，一般在200 ℃以内［14，18］，电阻率随着温度的升高而降低，电阻率与温度呈现负相关的变化关系。

图3 六块砂岩岩芯电阻率随温度变化（据参考文献 ［35］ 修改）Fig. 3 Variation of resistivity vs. temperature in six sandstone cores samples

直接利用曲线拟合出来的经验系数或者公式（2） 进行电阻率至温度的转换，来指导不同地质区地热资源的勘探，显然是不可行的。因为岩性、地层等因素的改变必然带来公式中经验系数的改变，使得特定的研究区域必须建立对小区块适用的经验系数。电阻率勘探方法及钻孔测温法在地热田勘查过程中是大量使用的方法，充分挖掘两者之间的内在联系，利用钻孔测温数据和电阻率数据对最终的温度预测结果形成有效的约束，是建立合理的经验系数的重点与难点。

国内相关学者基于可控源音频大地电磁测深的反演电阻率值，利用公式（2） 成功建立了有效的地下热储温度结构［16-17］。Huang 等人在此公式的基础上，充分结合钻孔测井电阻率、钻孔测温结果，以大地电磁测深的反演电阻率为基础，提出了通过深度进行分层约束的概念，并对电阻率进行归一化，建立了名为CCMOT 法（coefficient correction method of the optimal temperature，简称CCMOT）的电阻率和温度间转换的经验公式（3），在雄安新区地下热储温度结构预测模型的建立上取得了良好的效果［8］。随后其他学者利用CCMOT 法在雄安新区的其他大地电磁测量剖面上，进行了电阻率与温度间关系的转换，亦取得了较好的应用效果［19］。

式（3）中：T(x，z)—地下空间节点坐标A(x，z) 处的预测温度值，℃；RNT0(z)—不同深度岩石本身电阻率的归一化值，无量纲；αT0(z)—对应深度下的温度系数，无量纲；T0(z)—对应深度下的温度，℃；RNinv(x，z)—地下空间节点坐标A(x，z) 处电磁剖面上的电阻率归一化值，无量纲。

2.2 谷露地热田电阻率与温度的转换关系

借鉴经验公式（2）和（3）中的思路，通过电阻率比值和深度分层约束，建立了适用于谷露地热田的电阻率与温度转换公式，具体步骤如下：

1）根据音频大地电磁数据反演时剖分的不同深度的网格节点，对深度分层，并获得坐标A(x，Zi)。i即代表不同深度层的编号，Zi为第i个深度层对应的深度值。

2）选取BZ01 号钻孔（图2）作为此次经验公式建立时的标准孔，对钻孔的测温曲线进行线性回归分析（图4），得到深度大于1 100 m（标准钻孔最大深度）时温度估算的一元线性方程。即：

图4 BZ01 号钻孔测温曲线线性回归分析Fig. 4 Linear regression analysis of temperature measurement curve of borehole BZ01

式（4）中：TZi—正常背景下不同深度的温度估算值。

T(Zi)值的拟合，是为了获得无钻孔测温数据的深部地层在正常背景下的地温梯度。对整个钻孔的测温曲线进行拟合，因其包含热储层以及浅部盖层的高温数据，无法代表正常背景下的地温梯度，从而导致预测的深部地层的地温梯度误差变大。因此选取BZ01 钻孔测温曲线呈现典型线性变化特征的深度段进行曲线拟合。

（3）通过对公式（2）进行简单变形可得：

式（5）中：i—不同深度层的编号，无量纲；Zi—该深度层网格节点对应于地表的深度，m；αZ i—Zi处的温度系数，无量纲；Ts—给定的标准温度，℃；Rs—Ts温度下花岗岩本身电阻率，Ω·m；RZ i— 标准钻孔处AMT 反演电阻率，Ω·m；TZi—该深度所对应的钻孔测温值或估算值，℃。

需要注意的是，αZ i的计算中，如果对应深度有钻孔测温的实际值，则取实际值TZi，如果深度大于钻孔测温深度，温度值即取公式（4）中计算得到的TZi值。随后联合公式（2）、（4）和（5）得到：

式（6）中：R(x，zi)—地下不同位置在Zi深度处的AMT 反演电阻率，Ω·m；T(x，zi)—该位置对应的温度估算值，℃。

依据公式（6），即可由大地电磁反演得到的地下电阻率结构计算出同一剖面下的温度结构。

2.3 经验公式的验算及讨论

选取两条大地电磁剖面L1 线和L2 线（图2），其中L1 线穿过标准孔BZ01、验证孔YZ01，L2 线穿过验证孔YZ02。以标准孔BZ01 的钻孔测温和AMT 反演电阻率结果为基准，根据公式（6）计算验证钻孔处的温度值，并与实测钻孔测温曲线对比（图5）。

图5 实测温度与预测温度对比曲线Fig. 5 Comparison curve between measured temperature and predicted temperature

由图5 中的实测温度曲线与预测温度曲线可见，预测曲线的整体形态与实测曲线基本一致，且验证钻孔YZ01 的吻合度要优于YZ02。但在曲线中间段预测温度与实测温度均存在差别，分析其原因为温度值估算主要受电阻率比值这一因子影响，而AMT 反演电阻率无法完全等同于对应深度的岩石真实电阻率，仍然会受到上覆地层电阻率的影响，体现出一定的渐变性，因此当实测温度曲线在中间段波动较大时，温度预测值的误差变大。为了进一步对预测结果的准确性进行评价，引入拟合优度R2，从统计学的角度评价此次预测的准确程度。

式（7）中：TP—预测温度值，℃；TA—实测温度值，℃—实测温度平均值，℃。

R2的取值范围介于0～1之间，值越接近1，意味着预测值与实测值的拟合程度越好。经计算，图5a（YZ01）曲线与图5b（YZ02）曲线的R2值分别为0.823、0.710，表明此次温度预测值的精度相对较高，且验证钻孔YZ01 处的温度预测精度要高于YZ02，这亦与从曲线形态上的直观判断吻合。

3 剖面的电阻率结构

由电阻率结构导出温度结构，电阻率值是非常基础的数据。因此进行温度转换时，应先对电阻率的反演结果进行评价，以确保选用的电阻率值的有效性。

谷露地热田的基岩为高阻的花岗岩，低阻主要反映浅部第四系覆盖及断裂破碎带（图6）。在高程4 km 以浅，热田区主要的SN 向控热构造F4-1、F4-2和F4-3均能根据电阻率结构较为清晰的刻画，且F4-2断层有往深部下切花岗岩体的趋势，这一趋势在L1 线中反映得更为明显。这些特征与断裂带在地表的位置对应关系较好，且与该区之前的AMT 勘探成果［23-24，28］是一致的。因此选取的两条剖面的电阻率结果，可以用于进行温度换算。

图6 L1 和L2 线AMT 电阻率反演剖面Fig. 6 AMT reversely deducted resistivity profiles of L1 and L2 lines

4 剖面的温度结构及高温异常区预测

宏观上分析两条剖面的温度结构（图7），在高程4.3～3.9 km 这一深度带，均呈现西部高，东部低的特征；南北方向上来看，靠北的L1剖面上的温度值也大于靠南的L2 剖面。据此推断地下热流体在这一深度带有从北向南、自西向东的运移趋势。这一判断与谷露地热田实施的百米钻孔测温结果吻合［30］。

L1 剖面的预测温度结构中，在高程约4.3～4.0 km 范围内，可见明显的五处高温异常凸起，这五处异常，均与电阻率结果中的低电阻率异常一一对应，符合温度与电阻率呈现负相关性的特征；其亦与断裂构造的对应关系良好，符合该区深部热储主要沿断裂构造向上运移的主要特征［23-24］；YZ01 钻孔的热储层主要分布于高程4.3～4.05 km 范围，视厚度约250 m，而 Ⅲ 号高温异常凸起正对应于该深度范围，与实际的钻探成果吻合度很高；此次的温度预测结果与电阻率反演结果比较，温度预测结果也更为清晰地刻画出了这一热储层的位置及深度。据此推断，Ⅰ、Ⅱ、Ⅳ和Ⅴ号异常在高程4 km 以浅，存在热储的可能性较大，但其储层的厚度及最高温度均应小于Ⅲ号异常。高程4～3.5 km 范围分布的Ⅵ～Ⅸ号高温异常区，在深部热储的定位及探寻方向上，亦存在一定的指示意义：其中Ⅵ号异常已有钻孔BZ01揭露，深部温度高，空间分布的关系上与F4-1断裂更为紧密；Ⅶ号异常与F4-2断裂的空间分布密切相关，且有一定的由深部向上运移的趋势，推断该处温度异常存在较大可能为深部热储的反映；Ⅷ号异常深度较深，但温度较高，尚需结合地质资料进一步分析。

L2 剖面的温度结构中，高温异常位置亦对应于低电阻率异常，符合两者呈现负相关性变化的规律。在高程4.3～3.5 km 深度范围内，可见一明显的高温异常区 Ⅸ，空间分布上对应于F4-1、F4-2断裂，并具有往深部延伸的趋势。YZ02号钻孔揭露的热储层在高程4.35～4.13 km 范围，厚约220 m，位于 Ⅸ 号异常的边缘处120 ℃等温线的突变位置。L2 剖面对已知的温度异常反映较弱，整体对应情况不如L1 剖面，这从两个验证钻孔的拟合优度值上亦能说明。但从另一个方面来说，规模更大且温度更高的 Ⅸ号异常中心及深部，可能也具有更大的高温异常区，尚需进一步结合地质情况进行分析。

另外需要引起注意的是，虽然此次估算的剖面温度与已知的情况符合较好，并据此做出了有价值的预测，其主要基于谷露地热田地下岩性较为单一，主要为花岗岩，低电阻率与热储构造之间的相关性很强。但是目前仍然存在以下问题： 1）L2 剖面温度估算的结果差于L1 剖面，其原因可能是相关经验系数建立的标准钻孔BZ01 位于L1 剖面，离L2 剖面距离更远，不可避免地造成了地质因素及岩石物性因素差异上的增大，导致最终的预测结果精度下降，这意味着同样的经验系数值在空间的推广上存在较大的局限性。2）缺失了对应深度的钻孔测温数据约束，TZi也同时成为预测值，再来进行深部温度数据的估算，会进一步造成温度值误差的增大，尤其是在传热机制以对流为主的地热田中，甚至会出现深部温度值与常识相悖的现象。3）有效的电阻率数据的获得与标准钻孔的精细选取，对温度估算起到了决定性的作用，因此在不同的勘探区进行温度预测时，必须依据公式建立相应的经验系数，并应充分结合已知地质、钻探和地球物理成果，再来开展热异常区预测。

5 结 论

通过对已有地质资料的深入研究和分析，在结合前人研究成果的基础上，本文建立了适用于西藏谷露地热田的电阻率与温度转换经验公式，并进行了热异常区预测，得出以下结论：

1）电阻率比值及深度分层约束的理念，是此次建立电阻率与温度转换经验公式的基础。在以往研究中，电阻率比值主要采用钻孔测井电阻率作为比值项分子，而本文以钻孔处AMT 反演电阻率作为分子，具有一定的创新性和推广价值。

2）验证钻孔处的温度预测曲线与实际测温曲线形态基本一致，L1 和L2 两条剖面的温度结构具备由北向南、自西向东的运移趋势也与百米钻孔测温结果吻合［30］，基本证明此次温度预测结果的有效性。

3）在热异常区预测方面，3 号高温异常区与钻探揭露的热储层高度吻合，据此在两条剖面上圈定了四个浅部（Ⅰ、Ⅱ、Ⅳ和Ⅴ）高温异常区和三个（Ⅵ、Ⅶ和Ⅸ）深部高温异常区；其中Ⅶ 号异常与F4-2断裂的空间分布密切相关，且存在由深部向上运移的趋势，推断该处温度异常存在较大可能为深部热储的反映；而 Ⅵ号和 Ⅸ 号异常均与F4-1断裂密切相关，为下一步勘探提供了新目标。此外，本次温度估算所采用的相关经验系数值是在地下岩层主要为花岗岩的基础上建立的，理论上在不同勘探区域使用时，需要根据实际岩性和钻孔测温结果重新确定转换公式中所需的参数值。

致谢：感谢匿名审稿专家提出的建议。