粤港澳大湾区地热异常空间分布特征的GIS分析

徐讲湾， 黄少鹏， 魏正安， 李勇义

深圳大学土木与交通工程学院，深地科学与绿色能源研究院，广东 深圳 518060

能源问题已经成为人类社会经济发展中亟待解决的问题，开发利用可再生能源可减轻能源压力，并有利于抑制全球气候变化．地热能作为一种可再生能源，具有安全、稳定、利用系数高及综合价值高等特点［1］．据统计，包含粤港澳大湾区主体在内的广东省虽然是中国地热资源最丰富、温泉出露点最多的省份之一［2］，但目前广东省对地热能的开发利用力度尚显不足．缺乏对研究区地热资源成因机理和分布规律的认识是影响地热资源开发利用的重要原因之一．研究分析地热资源分布与地质构造和地球物理场之间的空间关联性，可以发现潜在的地热资源赋存区，降低勘探成本和风险，提高地热资源开发利用的经济效益．

通过地理信息系统（geographic information systems， GIS）集成各种地热地质、地球物理及地球化学等数据，并构建分析模型对地热资源进行空间分析的方法，已经在国内外地热勘探和资源评价中得到广泛应用［3-5］．现有的多种分析模型中，证据权重和加权信息量等数据驱动模型已被用于预测地质矿产资源分布［6］，且具有较好的预测效果［7］．目前，粤港澳大湾区地热资源开发利用的相关研究较少，利用GIS技术对粤港澳大湾区进行前期地热勘探分析的研究更是鲜有．

本研究系统收集和整理了粤港澳大湾区地热地质、温泉点分布、钻孔测温、大地热流、燕山期花岗岩空间分布、主要断裂、航磁异常、布格重力异常和居里面深度等数据，选取温泉及钻井地温正异常点作为地热资源点，构建证据权重（weight of evidence， WoE）模型和加权信息量（weighted information value， WIV）模型，定量分析了不同数据层与地热资源的空间关联性，并将不同数据层进行叠加分析来预测粤港澳大湾区地热资源的分布，以期为粤港澳大湾区地热资源勘探提供科学依据．

1 区域地质构造背景

粤港澳大湾区地层类型繁杂，从前寒武系到第四系均有出露（详情请扫描论文末页右下角二维码查看补充材料图S1）．地层序列由下至上可以概括为：① 深变质和强变形的前寒武纪结晶基底；②早古生代-中生代不连续的陆相碳酸盐岩单元；③晚中生代-第三纪大规模的岩浆侵入作用和火山活动冲破了的沉积层序；④ 第四纪海陆交互相沉积物［8-10］．自中生代以来，强烈的岩浆活动和构造活动在粤港澳大湾区内形成了一条宽400 km、呈NESW 向展布的中生代侵入岩带［10］，为粤港澳大湾区地热资源的形成提供了有利条件［11-14］．研究区侵入岩按其侵入时代主要可以划分为燕山期、印支期、海西期、和加里东期侵入岩，出露面积约占粤港澳大湾区总面积的三分之一，其中燕山期侵入岩面积最大．粤港澳大湾区从早白垩世末期至渐新世早期共有13期、123次以上火山喷发，主要分布于三水盆地［15］．区域中岩浆岩按岩性分类主要有中性岩和酸性岩等，绝大部分为酸性花岗岩．区域褶皱构造包括加里东期褶皱带以及华力西期-印支期褶皱带．根据地质及地球物理等资料分析，可以从西北往东南方向，依次把粤港澳大湾区划分为云开隆起区、阳春-花县盆地、增城-台山隆起区、斗门-惠阳盆地和粤东隆起区五个构造单元（图S1）．

粤港澳大湾区断裂构造的展布以NE-NEE、NW-NWW 和近EW 为主（图S1），主要发育于加里东期至燕山期挤压应力背景．NE-NEE向断裂带有：吴川-四会断裂带、郴县-怀集大断裂带、河源深断裂带、丽水-莲花山深断裂带、紫金-博罗大断裂和罗定大断裂带等；近EW向断裂带有：佛冈-丰良深断裂带、连平-恩平深断裂带和高要-惠来深断裂带等．其中，NE-NEE 向和近EW 断裂带为区域性控制构造，对区域侵入岩的空间分布起控制作用；NW-NWW 向断裂带形成时间相对较晚，且发育程度和出露规模相对较低，为次级断裂构造．NWNWW向断裂带往往与NE-NEE向和近EW断裂带相交，深大断裂与次级断裂的交汇为地表水的深循环提供良好的通道， 有利于热水、 热气的储存和运移．

2 数据选择和预处理

2.1 数据选择

通过对粤港澳大湾区温泉情况、钻井测温、大地热流、燕山期花岗岩分布、主要断裂分布、布格重力异常、航磁异常和居里面深度等数据进行综合分析，研究粤港澳大湾区地热资源的空间分布规律．

2.1.1 地热异常的认定

水热活动和高地温异常是地热资源的地表显示，统称为地热异常．因此，本研究以温泉点和高温钻孔作为存在的地热资源的标志，分析其与地热地质和地球物理场特征的空间关联．

① 温泉：目前发现粤港澳大湾区的温泉或地热田共41处．其中，14处温泉口温度为25～40 ℃；10 处为40～60 ℃；15 处为60～90 ℃；2 处为90～150 ℃．经统计，粤港澳大湾区温泉泉口平均温度约为50 ℃，为了细化地热异常程度与地热地质、地球物理场特征的空间关联性，本研究将泉口温度高于50 ℃的温泉点作为强地热异常点，低于50 ℃的作为一般地热异常点．

② 高温钻孔：通过钻孔测温可以圈定地温异常区，在热传导条件下，钻孔温度一般随深度线性增加；在某些浅部因素（如气候变化、土地利用变化和地下水活动等）的影响下，地温随深度的变化可能偏离线性［16］．对粤港澳大湾区部分钻孔进行测温，测点主要分布于惠州、广州和深圳等地，钻孔（如SZW-1 和NS-WQS-010）的温度-深度变化曲线见图1．

图1 粤港澳大湾区钻井测温曲线 （a） 所有钻孔； （b） 孔深0 ~ 100 m； （c） 孔深100 ~ 500 m； （d） 孔深500 m以上Fig.1 Borehole temperature profiles of the Guangdong-Hong Kong-Macao Greater Bay Area. (a) All boreholes, (b) borehole depth of 0 ~ 100 m, (c) borehole depth of 100 ~ 500 m, and (d) borehole depth greater than 500 m.

东南沿海地区的地温梯度为20～30 ℃/km［17］．根据广东省气象局、香港天文台和澳门地球物理暨气象局联合发布的《2021 年粤港澳大湾区气候监测公报》，粤港澳大湾区的2021 年平均地表温度为23.5 ℃，比常年偏高1.0 ℃．为了更直观地区分高温钻孔，本研究保守地以2021 年平均地表温度为起点，以30 ℃/km为斜率画直线（图1中虚线），作为判别高地温异常的依据，直线左边的曲线为正常地温或地温偏低的测温曲线，为非高温异常钻孔；直线右边的钻孔则为高温异常钻孔．根据图1 判别，粤港澳大湾区所测温的钻孔中高温钻孔共有16处（图2）．其中，惠热一井、ZK8、HSD-ZK001、JK-ZK1、JK-ZK4和YH-ZK005钻孔温度梯度不仅高且稳定，为强地热异常点，其余钻孔为一般地热异常点．

一般情况下，1 km2（1 km × 1 km）范围内的地热异常点属于同一地热系统［18］，选择一个点来代表该区域的地热异常点．本研究基于上述方法共识别出57 个地热异常点，其中，23 处强地热异常点，34处一般地热异常点（图2）．

2.1.2 大地热流

大地热流密度简称大地热流，指单位时间、面积内，地球内部以热传导方式传输到地表的热流．大地热流可以综合反映一个地区的地热场特征，是了解地热异常区分布规律的重要数据．本研究的大地热流数据来自国际热流委员会（International Heat Flow Commission）热流数据库中的中国及邻区的大地热流测量值和姜光政等［19］的中国大陆地区大地热流数据．经过筛选处理及反距离加权法插值，粤港澳大湾区及邻区大地热流图如图3（a）．

图3 本研究选用的数据层 （a）大地热流；（b）燕山期花岗岩分布；（c）距断裂距离；（d）居里面深度Fig.3 Data layers selected for this study. (a) Heat flow, (b) Yanshanian granite, (c) distance to a major fault, and (d) Curie depth.

由图3（a）可见，粤港澳大湾区大地热流值为60.3～88.9 mW/m2，平均值为72.4 mW/m2，高于中国大陆地区平均大地热流值61.5 mW/m2［19］．在肇庆西部、江门西部、惠州、深圳及香港等区域表现出大地热流值较高的特征，以及由北向南热流升高的趋势．

2.1.3 燕山期花岗岩

粤港澳大湾区燕山期中酸性岩浆侵入活动强烈．花岗岩的放射性生热率高，其广泛出露可能与地热资源的空间分布密切相关．本研究将粤港澳大湾区构造图进行矢量化，从而得到燕山期花岗岩的分布，如图3（b）．粤港澳大湾区广泛存在的燕山期花岗岩不仅是潜在的热源体，且花岗岩体边部因冷凝收缩可为地热资源赋存提供有利的条件．通过圈定燕山期花岗岩及附近一定范围区域，有利于寻找隐伏地热异常．运用欧氏距离分析计算粤港澳大湾区中每个点到花岗岩的距离，再通过空间叠加，分析地热异常点与燕山期花岗岩分布的空间关联．统计发现，约65%的地热异常点处于燕山期花岗岩体内或边界上，约93%的地热异常点处在距离燕山期花岗岩体边界3 km范围内．

2.1.4 断裂

断裂活动与地热资源分布联系密切，断裂活动形成的裂隙为地下热水深循环提供了通道［20］．粤港澳大湾区断裂构造发育，展布方向主要有NE-NEE、NW-NWW 和近EW 向．总体上，压扭性断裂对地热水的储存及运移不利，但构造复合部位，可能具有张性或张扭性的破碎带，因此温泉大多出露在断裂构造的复合部位，在距断裂较近的区域，存在地热异常现象的概率较大．通过对断裂构造分布进行欧式距离分析，可得到区域每个点到断裂距离，如图3（c），从而来量化地热异常点与断裂空间分布的关系．在距离分析结果中，每个栅格属性值代表与其最近断裂的距离．统计显示，约76%的地热异常点分布在主要断裂的10 km范围内．

2.1.5 布格重力异常

研究区布格重力异常由北向南呈正负交替分布，四会-花都-增城-博罗以北为布格重力负异常区，负异常区强度变化幅度在35 mGal左右；四会-花都-增城-博罗与开平-新会-中山-光明之间为正异常区，该正异常区跨越佛山、广州黄埔、东莞和惠州惠城，重力异常值最高为16 mGal．开平-新会-中山-光明以南为布格重力负异常区，该负异常区强度变化幅度小，变化幅度在10 mGal左右［21］．

引起布格重力异常的地质因素主要有侵入岩和莫霍面埋深［22］．重力场分解通常使用各种频率滤波方法．XI 等［23］对广东省布格重力异常进行了小波多尺度分解，利用不同的细节阶数和四阶近似来估算不同布格重力异常源的埋深．分析结果显示，广州及其周围地区的高异常来源于地壳深部，与莫霍面隆起有关；区域内低重力异常与花岗岩有明显的关系．任镇寰等［21］对珠江三角洲地区的布格重力资料进行了上延计算分析，描述了区域内深大断裂的结构特征．结果显示，珠江三角洲的布格重力异常梯度带与区域的断裂带分布密切相关．布格重力异常梯度高，表明地下可能存在有利于地热异常的断裂构造，存在形成地热异常的有利条件．因此，通过布格重力异常可以圈定区域性断裂构造及与地热系统相关的局部热异常．

2.1.6 航磁异常

磁异常分布常与地下水热活动及区域构造应力变化有关．在热储区，高温会使围岩发生热蚀变作用，导致周围岩体的磁性减弱或产生退磁作用［24］；另外，构造应力会使得岩石的磁性在应力方向减弱．航磁异常低值区域较高值区域存在地热资源的可能性大．粤港澳大湾区航磁异常值变化区间为-35～52 nT，最大和最小磁场值分别位于广州从化和惠州博罗地区．粤港澳大湾区的磁场具有明显的分区性，以鹤山-顺德-南沙-博罗一线为界．界线北部磁场较南部变化平缓，以负磁场为背景，局部分布有NE、NW 向正磁场．界线南部磁场以正磁场为背景，正负磁场交替变化为特征［25］．把航磁异常分布和地热异常点进行空间连接，结果显示，约87%的地热异常点分布于航磁异常值小于5 nT的区域内．

2.1.7 居里面深度

居里等温面简称居里面，是地球岩石圈上部的一个温度界面，指岩石圈中的磁性岩石矿物由于高温而失去磁性的温度界面．居里面深度与大地热流、地温梯度、断裂活动及岩石物性等因素有关［26］，对认识区域深部热状态有重要意义．一般来说，较浅的居里面深度对应较高的平均地温梯度，居里面深度的分布与地壳厚度、岩浆岩侵入活动有关［27］．因此，可以将居里面深度分布作为研究区域背景地温场的重要参考数据．粤港澳大湾区及邻区居里面深度分布如图3（d），数据由原广东省地矿局物探大队根据全国统一规范实测和编制的广东省1∶5 × 105航磁异常图进行矢量化和计算得到．粤港澳大湾区居里面埋深在12～22 km，佛山中部、广州从化区及深圳等地存在明显隆起现象．中国陆域的居里面深度为18～40 km，平均埋深约为30 km［27］，粤港澳大湾区居里面整体埋深较浅．

2.2 数据重分类

对数据层进行重分类是进行空间数据叠加分析的重要步骤之一，根据各数据层的阈值进行合理的分级，划分的类别越多越有助于提升结果的可靠度［3］．目前，对数据层进行重分类的方法有等间隔分类、自然间断法分类等［4-5，7，28］，本研究采用自然间断法将各数据层进行重分类，请扫描论文末页右下角二维码查看补充材料表S1重分类结果．

另外，还需要对数据层栅格单元大小进行网格划分，网格划分的越小，每个单元中的地质环境越相似，越有利于提高模型分析的精准度．汤国安等［29］考虑研究区面积、尺寸比例及计算机性能等因素，提出了网格大小划分标准．粤港澳大湾区的面积约为5.6 × 104km2，地图比例尺为1∶15 × 105，将粤港澳大湾区数据层的网格大小确定为100 m ×100 m，所有数据层被划分为5 509 980 个网格单元．重分类和网格划分的结果可扫描论文末页右下角二维码查看补充材料图S2．

3 数据整合模型构建

将多维空间数据进行整合的方法主要有知识驱动模型和数据驱动模型．它们的区别在于：知识驱动模型的参数为专家经验和评估得到，例如布尔逻辑、指数叠加和模糊逻辑模型等；数据驱动模型的模型参数是依据具体的数据进行计算得到，如证据权重、加权信息量模型和逻辑回归模型等［30］．粤港澳大湾区地热资源成因复杂，专家经验评估得到的参数难以表征研究区地热资源空间分布特征，而区域内的温泉和高温钻孔等具体数据则可以用于计算模型参数．因此，选取数据驱动型模型对粤港澳大湾区数据层进行整合．

3.1 WoE

WoE是一种基于贝叶斯概率框架的空间分析方法，可以量化数据层与地热异常点之间的空间分布关系，从而分析数据层与地热异常分布规律的关联性，确定不同数据层的不同分类对地热异常的权重［31-32］．获得权重后，再对不同的数据层进行空间叠加，研究地热异常空间分布规律．构建WoE 模型的步骤如下［31］．

不同的数据层与地热异常点的空间关联性可以用正、负权重（W+和W-）表示，

其中，Np和分别表示在研究区中不同数据层存在地热异常和不存在地热异常的单元个数；G和-G分别表示在研究区中存在地热异常点和不存在地热异常点的个数．

W+> 0和W-< 0均表示数据层的某个分类与地热异常点的空间关联性呈现正相关，在该情况下，出现地热异常的概率较大；反之，该数据层分类与地热异常点呈现负相关，即出现地热异常的概率较小．数据层i的某个分类j与地热异常点之间的空间相关程度可以用Cij表示，Cij> 0 为正相关，Cij= 0为不相关，Cij< 0为负相关．

为减小偶然因素对Cij的影响，可以对Cij进行标准化对比度S（C）计算，S2(W+)和S2(W-)分别为W+和W-的方差，为

其中，N{}· 为数据层与地热异常点叠加时的点数．通过选取Cij/S(C)的最大值来确定各数据层的最优分类阈值，并根据最优阈值对数据层中的单元进行二值化处理，数据层单元格属性值在最优阈值内赋W+值，其余赋0值．最后，将二值化后的数据层按式（8）进行叠加处理，得到地热异常分布的预测图．

其中，P(G)为先验概率；n为选取的总数据层数；N为研究区总单元格数；V为栅格单元输出值．

对地热异常分布规律进行空间叠加预测分析时，如果各证据层之间存在较强的相关性，则会重复考虑数据层与地热异常空间分布的关联，从而导致模型结果与实际情况产生较大误差．为减少该误差，需要对数据层进行条件独立检验，本研究通过卡方（χ2）检验来判断不同数据层之间的相关性［30］．经计算，各数据层之间的相关系数均小于阈值3.84，请扫描论文末页右下角二维码查看表S2．虽然航磁异常与燕山期花岗岩的χ2值接近阈值，表明它们存在一定的条件依赖性，但是在5%的显著性水平的假设下认为满足条件独立，可以作为模型分析计算的数据层．

3.2 WIV

地热异常的出现经常伴随多种地热地质、地球物理场的异常出现，在多种异常特征下总会存在某一种最优组合，使地热异常出现的概率最大［3］．WIV则可以量化地热异常与地热地质、地球物理数据层的关联密切程度［33］．通过把与地热异常有关的特征或现象转化为信息量值，信息量值越大，发现地热异常的概率越大，反之，地热异常出现的可能性较低．构建WIV的步骤［33］如下．

不同数据层i的信息量值用Ii表示，将重分类后的数据层按式（10）和式（11）计算可得到不同数据层的信息量值，

其中，Nij为第i个数据层第j分级中地热异常点数；N为总地热异常点数；Sij为第i个数据层第j分级中栅格单元数；S为总栅格单元数；m为第i个数据层的分类数．

不同数据层与地热异常的空间关联可能存在差异性，通过对数据层进行加权处理，可以评估数据层的差异性对模型分析结果的影响．权重系数Wi可由式（12）至式（15）［34］确定，

其中，Gij为地热异常点密度；Pij为第i个数据层映射的第j分级的归一化值；Hi为信息熵；n为数据层数．栅格单元的加权信息量值I为

4 模型分析结果

4.1 数据层空间关联性

为表征不同数据层与地热异常的空间关联性，将粤港澳大湾区地热地质、地球物理数据与地热异常点的空间分布特征进行量化，量化结果可扫描论文末页右下角二维码查看补充材料表S3．证据权重模型考虑了同一个数据层的不同分级与地热异常现象的关联性．空间关联性分析表明，大地热流值大于75 mW/m2、距离断裂小于3 km、布格重力异常梯度大于1.6 mGal/km、燕山期花岗岩边界或内部、航磁异常值为-5～5 nT及居里面深度小于15 km的区域与地热异常的空间分布较为密切，粤港澳大湾区的大部分地热异常点都出现在上述区域中．

由于不同数据层与地热异常空间分布的关联性存在差异性，在进行空间叠加分析时，需考虑不同数据层之间的差异性以提高结果的准确性．而加权信息量模型正是考虑了这一点，量化结果显示居里面深度、燕山期花岗岩和主要断裂数据层与粤港澳大湾区地热异常空间分布的关联性较强，航磁异常数据层与粤港澳大湾区地热异常空间分布的关联性不强．因此，在圈定地热异常区的时候应重点考虑居里面深度、燕山期花岗岩和主要断裂这几个数据层．

4.2 地热异常程度与数据层空间关联性

在论述地热异常程度时，将地热异常点划分为强地热异常点和一般地热异常点．如果地热地质、地球物理等特征与地热异常的空间关联性较强，则强地热异常现象发生的概率较大．

为进一步验证4.1节中与粤港澳大湾区地热异常空间分布关联性较强的不同数据层的最优组合，把4.1节计算得到的各数据层与地热异常的空间关联性量化值用自然间断分类方法分为低、中、高空间关联性3类，结果请扫描论文末页右下角二维码查看补充材料表S4，然后分别统计不同程度的空间关联性区域内强地热异常发生的概率．为减少偶然因素对统计结果的影响，将强地热异常点的数量除以每类区域中栅格单元数与该区域内总异常点数的乘积．如图4，在不同数据层的高空间关联性区域，强地热异常点的占比均最大， 换言之，大地热流值大于75 mW/m2、距离断裂小于3 km、布格重力异常梯度大于1.6 mGal/km、燕山期花岗岩边部、航磁异常值为-5～5 nT 及居里面深度小于15 km 等区域，可作为与粤港澳大湾区地热异常空间分布关联性的最优组合．因此，选择上述组合进行空间叠加分析以确定粤港澳大湾区地热资源的空间分布特征．

图4 不同关联程度区域中强地热异常点概率Fig. 4 Proportion of strong geothermal anomalies in different spatial relevance categories with low spatial relevance region （blue） and medium spatial relevance region （gray） and high spatial relevance region （red）.

4.3 模型结果重分类

证据权重模型和加权信息量模型的计算结果如表S3．为更直观地揭示粤港澳大湾区地热异常空间分布特征，需要对模型按计算结果进行重分类［3］．采用自然间断法将证据权重模型和加权信息量模型的计算结果分为5类．在加权信息量模型计算结果中，信息值越小于0表示地热异常现象不存在的可能性越高，信息值等于0表示地热异常现象存在的概率几乎为0，信息值越大于0 表明地热异常现象存在的概率越高．据此，将信息值0 设置为第1 个阈值，将加权信息量模型的计算结果小于0的部分分为1类，将大于0的部分按自然间断法分成4类．基于证据权重模型和加权信息量模型的重分类数据层空间关联性分析表明，粤港澳大湾区内高潜力地热异常区的位置在空间上基本重叠，主要分布在江门的西部，肇庆中部，广州、惠州的北东以及深圳等地（图5）．高潜力地热异常区主要位于低山或丘陵区，燕山期花岗岩分布广泛，且热导率较高，易出现强地热异常．相反，低潜力地热异常区多位于沉积盆地，上覆低热传导率的沉积盖层（主要为第四系海陆交互相沉积物），限制了热量传输，不易形成强地热异常．以上分析与模型结果也是相吻合的．

5 结 论

基于GIS 技术，选取粤港澳大湾区地热地质、地球物理资料，通过建立证据权重和加权信息量模型，对大湾区地热异常区进行叠加分析，得出以下主要结论：

1）居里面深度、燕山期花岗岩中高含量的U/Th/K放射性元素生热率和主要断裂数据层与粤港澳大湾区地热异常空间分布的关联性较强，航磁异常数据层与粤港澳大湾区地热异常空间分布的关联性较差．

2）与地热异常有关的地热地质、地球物理场的数据层组合中，最优组合为大地热流值大于75 mW/m2、距离断裂小于3 km、布格重力异常梯度（即变化率）大于1.6 mGal/km、燕山期花岗岩边界、航磁异常值为-5～5 nT及居里面深度小于15 km的区域．

3）粤港澳大湾区存在隐伏地热资源的可能性较高的地区包括江门西部、肇庆中部、广州东北部、惠州以及深圳等地．