0 引言

青藏高原东南缘位于欧亚和印度两大板块会聚、消减、相互作用的边缘地带，是中国大陆地壳运动最强烈、地震活动最频繁的地区之一，这种特定的地质环境是局部高温地热带形成的基本地质条件(陈墨香等，1994).地层深部物理特征变化与深部地热活动有着密切关系，地热活动或热状态的改变对孕震区介质的物理和化学性质都有着重要的影响，是触发大地震的重要因素(王云等，2018).因此，地热学分析，不仅是开发利用西南地热资源的重要基础，也是分析该地区构造活动背景与地震活动等的有效途径.

大地热流，指由地球内部传输至地表并释放至空中的热量，它是表征地球热状态的重要地球物理参数，是地下热能在地壳表层的综合反映，被诸多学者称为“窥视”地球内热的窗口(郝春艳等，2014).获取热流数据最直接有效的方法是利用钻井技术，通过地温测井技术测得地下地层的地温梯度，经过取心、采样，在实验室内测得岩石热导率，进而计算出该深度下的热流值.然而，高质量的热流测量在全球范围内是很难达到的，地表热流也易受到局部热异常的污染(Ross et al., 2006)，如断裂分布、地下水循环、地形变化和气候变化等，精确的热流数值需要深部钻探获取，难度较大.青藏高原东南缘的地表热流数据稀少，实测数据仅覆盖小部分区域(如图 1)，Huang和Zhao(2006)计算了中国大陆地区深部P波扰动的径向剖面，其径向剖面揭示出腾冲地区的低速异常延伸至深部200 km左右，可推知该地区地热异常的影响主要来自深部；王云等(2018)发现滇东红河地区地壳深部20~40 km范围内广泛存在高导低速层，推断可能富含高温高压流体或部分熔融.由于缺少地壳内部热流数据的约束，无法直接确定地球内部的地热状态，阻碍了深部动力学及构造运动的研究.

居里面是指磁性矿物在一定温度下从铁磁性转变为顺磁性的等温界面，是磁性层的下界面.居里面作为岩石圈内的温度界面，与地表热流的分布有着非常密切的关系，根据一个地区的居里面深度变化结合生热率和热导率等参数可以得到地层热流分布等有价值信息(Elbarbary et al., 2018)，因而可用于估算青藏高原东南缘热流分布.在前人的诸多研究中，根据居里面深度求取大地热流的方法大致分为两类：一类是假设地层为无热源层，地层热导率均匀，遵循简单的一维傅里叶热传导定律(Nwankwo and Shehu, 2015；Saibi et al., 2015；Li and Wang, 2016；Saada，2016；Afshar et al., 2016；Bello et al., 2017)；另一类是假设地壳中存在热源(即放射性元素衰变产热)，地层热导率均匀或纵向分层变化(Lowry and Pérez-Gussinyé，2011；汪集旸等，2012；Speranza et al., 2016；Martos et al., 2017).第一类方法相对来说计算粗略，第二类方法比第一类方法考虑因素更全面，更接近真实热流值.

本文在第二类方法的基础上，综合考虑深部热源向上传导和地壳放射性元素衰变产热的影响.首先，利用地热参数之间的关系分别迭代计算出研究区地壳上下层的热导率分布，突出了热导率对大地热流数值的影响，然后，从居里面深度出发，计算了青藏高原东南缘的大地热流分布，最后得出地下不同深度处的热流分布，并归纳出其与地震活动之间的联系.

1 区域地热背景

青藏高原东南缘地处地中海—喜马拉雅地热活动带最东端，属于陆-陆碰撞型地热带，该地区是几大地质构造单元汇聚的地区，这些块体的发育历史虽不尽相同，但多期性、持续性的地质构造活动从未停止，特别是新近纪地史时期构造活动表现剧烈，活断层高度发育，地震频度高、震级大而分布很广.这种活动性程度高的地质条件也是该地区具有较高地热背景值的基本原因.该区域整体上呈现沿北东方向逐渐降低的地热背景，这种大地热流变化趋势反映了青藏高原“东构造结”北东向推挤的热传递效果(陈墨香等，1994；张健等，2017).

青藏高原东南缘的局部地热异常分布也十分明显，川西地区多处出露温泉，如康定、理塘、乡城和甘孜等地，这些水热活动区伴随断裂带的发育而展布，分别沿着鲜水河断裂带、甘孜—理塘断裂带和金沙江断裂带(如图 1).滇东地区的红河为典型地热田，该地区的红河断裂带为活动性较强的深大断裂，为浅层地热发育提供了条件.滇西地区的腾冲火山是我国仅有的几个活火山之一，该区地温梯度高，高温热泉密度大，数目多，而且类型复杂齐全，有高温沸泉、温泉和地热蒸汽等(徐世光和郭远生，2009).

通常地表热流分布还受地壳内部基底起伏与局部熔融的影响，这些干扰因素来自地球内部，往往与幔源热流体的向上运移有关.此外，地热与地震活动之间存在关联.地震活动在地热的形成过程中起了一定的作用，反之地下热通量的形成和运动过程中所引起的一系列物理、化学变化，对于地震的孕育和发动也给予比较大的影响，地震是热异常的表现(李杨等，2018)，地热是地震在应力上的体现.越来越多的学者认为大多数地震事件不发生在热流高值点，而是发生在高低值过渡带上，尤其是大震级事件.

2 数据 2.1 居里面深度数据来源

本文数据来源于杨海和熊盛青等人计算的中国陆域居里面深度图(杨海，2015)，该数据是基于最新汇编的1:100万中国陆域航磁ΔT化极数据，采用功率谱法计算而来.将其与侯重初版的全国居里点深度对比，显示两者居里面特征一致，其中康滇地区的居里面深度与申宁华等(1986)居里点数据对比显示两者趋势一致，深度较一致，相比前人的计算结果，该数据的分辨率更高、细节更多.此外，通过与实测热流数据对比发现，高热流值区对应居里面深度隆起地区，低热流值区对应居里面深度凹陷区(如图 1)，符合一般性地质与地热学规律.本文从中矢量化出青藏高原东南缘(22°N—35°N，95°E—105°E)的居里面深度数据(如图 2)，横向网格间距为0.15°.该地区居里面起伏明显，深度在18~44 km之间，隆起区呈南北向，分别沿玛沁—德格—乡城—昆明—红河一带和玛沁—德格—乡城—腾冲—临沧一带展布，凹陷区以玛沁—甘德的西部与东部两侧、四川德阳与云南楚雄及以南地区为主.

2.2 放射性生热元素数据来源及地表产热量计算

地表热流与放射性生热元素的丰度有着很强的相关性，大陆大约一半的表面热流可以归因于大陆地壳中铀、钍和钾的放射性同位素的产热.在自然界中，最主要的放射性核素有铀(²³⁸U)、钍(²³²Th)和天然钾中不稳定同位素⁴⁰K.其中，铀元素有²³⁵U、²³⁸U和²³⁴U三种同位素，²³⁴U和²³⁵U的含量远远小于²³⁸U，三种同位素共生在一起，因此，可以将铀元素的丰度值视作²³⁸U核素的丰度值.钍元素的同位素较多，但是天然钍元素几乎全由²³²Th组成，因此，钍元素的丰度值可以视作²³²Th核素的丰度值.本文使用的铀钍元素丰度数据来自中国地质科学院地球物理地球化学勘查研究所于2012年发布的全国铀地球化学图和全国钍地球化学图.钾元素有三个同位素，但只有⁴⁰K有放射性，所以，本文将铀(U)和钍(Th)的丰度值C_U和C_Th分别赋予一半的权重，利用以下关系计算出钾(⁴⁰K)的丰度值C_k(Turcotte and Schubert, 2002)：

(1)

青藏高原东南缘地区的铀(²³⁸U)、钍(²³²Th)和钾(⁴⁰K)的丰度值分别为图 3a、b和c.由此，地表产热量A可根据以下经验关系式计算得出(汪洋，2006)：

(2)

其中，H^U=9.52×10^-5，H^Th=2.56×10^-5，H^K=3.48×10^-5，分别代表铀钍钾元素丰度值的系数.H₀指地表单位质量岩体的生热率，单位W·kg^-1.密度ρ为近地表平均密度，本文选取密度值为2.3 g·cm^-3(唐新功等，2012).

2.3 实测大地热流数据

实测热流数据来自中国大陆地区大地热流数据汇编第二版、第三版和第四版(汪集旸和黄少鹏，1990；胡圣标等，2001；姜光政等，2016)，其中，研究区域内的热流数据共77个，质量级别为A、B类的数据共58个，本文选取此58个数据，用来约束地热相关参数的计算，以及与估算热流值作对比.

3 方法原理 3.1 热流估算方法

本文采用一维稳态热传导地层模型，认为地球深部的热源和地壳放射性生热元素产热向地表作径向传导.在稳定的大陆地区，表面热流随着地表岩石的年龄而系统地减少，同样地，岩石中放射性同位素的浓度也随着岩石年龄的增加而减少，但是只有在时间尺度较大时才会引起明显的热流变化，因此，本文不考虑温度和热流在时间上的变化，并且不计入浅层的地下水、断层或裂隙发育带的渗透等带来的热影响.

实际地层中放射性核素产热不均匀，通常认为其产热量A随深度z的增加呈指数衰减，并且与地层温度T存在以下关系(Lowry and Pérez-Gussinyé，2011)：

(3)

其中，k代表热导率，单位W/(m·K)，T代表地层温度，单位℃，h_r指放射性生热元素富集层的厚度，单位km(详细说明见后文).对深度z作一次积分，并代入边界条件(Turcotte and Schubert, 2002)：

(4)

此处q_m为地幔热流，单位mW·m^-2.得到：

(5)

通常来说，居里面温度T_c为580 ℃，近地表z₀温度T₀为15 ℃，根据(5)式对深度作二次积分，推算出温度T(z)为

(6)

对(6)式作进一步变换即可得到大地热流q₀与居里深度z_b的关系(Martos et al., 2017)：

(7)

利用大地热流q₀可以计算地下其他深度z处的热流q(z)分布(Turcotte and Schubert, 2002)：

(8)

当需要根据岩性或地层划分而赋予热导率不同数值时，地表热流和地下不同深度处温度的计算会更复杂一些，本文在Andrés等(2018)的基础上分别对公式(6)、(7)作出改进，得到当热导率分为两层时对应模型下的地表热流和地下不同深度处温度的计算公式，如下：

(9)

(10)

其中，z代表当前计算深度，k₁、k₁分别代表上下地层的热导率，h₁、h₂分别代表上下地层的厚度.

3.2 放射性生热元素富集层的厚度计算

放射性生热元素富集层的厚度是指在该深度处地层放射性产热量是地表处的1/e.为了确定研究区域的放射性生热元素富集层的厚度，需要结合实测热流数据和地表产热量，根据两者关系作线性拟合，其直线的斜率即为放射性生热元素富集层的厚度h_r.具体关系为(Lowry and Pérez-Gussinyé，2011)：

(11)

由于研究地区的实测热流数据较少，为提高参数h_r的准确度，本文在公式(11)的基础上作了多次反复拟合，原理如下：利用实测热流值q₀和地表产热量A作线性拟合可得到斜率值h_r¹，再由斜率值h_r¹反推得到新的热流数据q₀¹，不断重复此过程，直到qⁱ₀与q₀^i-1之差满足误差要求，确定最终厚度值h_r.

3.3 热导率计算

岩石热导率表示地壳岩石的导热能力性能，地层中温度、压力、各向异性体和孔隙性等都会对岩石的热导率值产生影响.研究地区的居里面深度起伏较大，这导致居里面以上地层的地温梯度在横向上变化明显，采用热导率平均值的做法势必会压制热导率参数对局部热流值的影响，因此，本文将热导率分为上下两层，并考虑了温度变化对热导率数值的横向影响，忽略了对热导率影响次之的地层压力因素.

假设地表到居里面之间的地层为一系列大小相等的长方体组合，每一长方体内部温度均匀，热导率随温度变化的关系为(Lowry et al., 2000)：

(12)

其中，a、b值参考了Lowry的经验值，a=0.28；b=3.16×10^-4.利用地热参数之间的关系作多次迭代计算，可以得到研究地区热导率的分层分布.以研究地区的某一网格点为例，具体计算流程及流程图如图 4.

(1) 选取适当的热导率初始值，代入单层热导率模型下的地表热流计算公式(7)，得到初始的地表热流估算值.并将上述地表热流估算值代入地层温度计算公式(6)，得到地下不同深度处的温度分布.

(2) 根据公式(12)，在获取地层温度的基础上可以得到该地层的热导率分布.并根据上下地层的深度分界线，求出上下地层对应的热导率平均值k^up和k^down.

(3) 将得到的热导率平均值k^up和k^down代入双层热导率模型下的地表热流计算公式(9)，然后计算地下不同深度处的温度分布，如果计算深度小于热导率分层深度h₁，则将得到的热导率平均值k^up代入单层热导率模型下的地层温度计算公式(6)，如果计算深度大于热导率分层深度h₁，则将得到的热导率平均值k^up和k^down代入双层热导率模型下的地层温度计算公式(10).

(4) 重复步骤(2)—(3)，直至上下地层对应的热导率平均值k^up和k^down与前一次计算结果之差小于误差允许范围，结束计算，并得到最终的地表热流估算值.

本文的具体参数设置为：热导率初始值为2.51 W/(m·K)(汪缉安等，1990)，每个长方体厚度为1 km，热导率上下层分界面为地下15 km(详见下文4.1节)，每次热导率迭代的允许误差为小于0.01 W/(m·K).

4 结果与讨论 4.1 研究区放射性生热元素富集层厚度及热导率分布

从整个研究区域来说，放射性生热元素富集层的厚度h_r变化的范围不大，本文最终拟合出研究地区的放射性生热元素富集层厚度h_r为8.9007 km(见图 5).在大多数情况下，放射性生热元素富集层的厚度h_r在10 km左右，且对热表热流的数值计算影响不大.

本文根据研究区的地球物理探测结果与物性分布特征分层设定了热导率参数.横跨滇中地块的大地电磁探测剖面显示腾冲—楚雄一带地壳的电性结构主要分为两层，上层为高阻，下层为低阻，分界线在地下15 km左右(于常青等，2017)；四川西部松潘—德阳一带的大地电磁测深结果显示该地区地下10~15 km处出现明显的电性转换(朱迎堂等，2008)；

云南地区的背景噪声层析成像结果显示该地区地下15 km处为S波速度变化分界面，界面以上S波波速较低，界面以下S波波速较高(秦维秉，2017).因此，本文将地下15 km设定为青藏高原东南缘参数热导率分界面，最终计算并得出青藏高原东南缘的热导率分布(如图 6)，上层热导率平均值为3.04 W/(m·K)，下层热导率平均值为2.31 W/(m·K).该地区的西南边缘及康定—九龙等地为热导率低值区，这与该区域的居里面深度较浅、地表产热量较高的特点相对应，符合温度较高时热导率参数较低的规律.

4.2 研究区热流分布

在获取研究地区的热导率分布和放射性生热元素富集层的厚度数据的基础上，本文估算出青藏高原东南缘的地表热流分布(如图 7a)，青藏高原东南缘的大地热流数值位于44~108 mW·m^-2之间，平均值为74 mW·m^-2，符合研究地区西南高、东北低的背景趋势.根据大地热流的分布进而得到了地下10 km、20 km和30 km深度处的热流分布(如图 7b、c和d).随着深度的增加，研究地区内地壳内部的热流值整体降低，幅值区间缩小，增幅趋于缓和.其中，康定及其周围地区的热流值降幅明显，说明该地区的热源并非来自深部，而是浅层的局部热影响.这里需要强调的是，热流值的大小不等同于温度的高低，影响热流值的因素有岩石热导率和热源的大小，热流低值可能是由于热源的减少，或者是由于温度随深度增加使得岩石热导率降低，综合该地区的地热背景，本文认为是两种因素的共同影响.

本文将地表热流估算值与实测数据进行对比，以分析结果的可靠性.图 8显示，超过一半的数据点与实际数据误差在10 mW·m^-2以内.其中误差较大的散点均体现出实测热流数据小于本文估算值的特点，这是由于实际测井的深度在地下100~7000 m之间不等，而本文的结果指示地表值，因此实测值小于等于估算值是合理的.此外，在实际地层中，地表热流受断裂分布、地下水循环、地形与气候等的影响，估算热流值未计入此类热干扰因素的影响，导致二者存在一定差异.其中，基底起伏与局部受热熔融产生的热影响是包含在本次估算热流值之内的，当地下基底受热隆起或地壳内局部受热熔融时，地层温度会发生变化，对应地区的居里面起伏状况也会随之改变，进而将信息携带在居里面深度数据中.因此，地表热流的估算值包含的地壳浅部信息较少，包含地壳深部的信息较多，更有利于理解岩石圈深部热结构.同时，相比密度、地震波速等其他地球物理参数，大地热流数值的不稳定性较大，综合以上，虽然实测数据具有一定的可靠性和准确性，考虑到浅部干扰因素大和深度的不一致性，实测数据对于本文估算值而言不具有绝对否决性.

为了进一步验证数据的可靠性和合理性，本文按照数据散点的区域分布将其划分为点簇，以分析数据在空间分布上的差异.对比可知，实测热流数据(图 9a)在相邻测点之间数值具有跳跃性，估算热流值(图 9b)的数据连续性更好，并且，后者体现出了数据的变化趋势：整体上由南向北降低，在北纬29°南侧地区，由西向东先降低再升高.因此，认为本文所采用的数据及参数设置合理，得出的大地热流值可靠，对青藏高原东南缘地区适用.

相比大地热流的数值大小，热流的异常高低与所在的区域的分析更具实际意义.通常，大陆地区的平均大地热流值在60 mW·m^-2以下，超过80~100 mW·m^-2即显示出明显的热流异常.本文将热流的高低异常分布与前人的研究结果对比，发现两者在大多数地区具有一致性(见表 1)：在四川西部，热流高值分布在金沙江断裂、德格—乡城断裂、甘孜—理塘断裂沿线和康定等地区，沿西北-东南方向展布.在云南地区，热流高值区分布广泛，分为滇东和滇西两区：滇西地区受地质构造的隆起、拗陷布局的控制，以腾冲地热区为典型代表，分布有临沧等热流高值区；滇东地区的热流高值区主要分布在滇东断陷盆地，如红河、昆明以及小江断裂带沿线.在四川东部与云南的楚雄等地为热流低值区.通过对比，与前人研究结果不一致的为康定北段和小江断裂以东地区，可解释为：康定北段地区的浅层断裂促进了热流的横向传导，进而扩大了热流高值区的分布；小江断裂以东的热流高值范围分布较广，同时，该地区居里面深度较浅，地壳浅层产热量较高，这些特点均与地表高热流值相符合，考虑到该地区的实测数据虽为低值但点数稀少，因此，认为小江断裂以东为热流高值区.

4.3 地壳内热流分布及与地震活动的相关性

本文将研究地区不同深度处的热流值与青藏高原东南缘的地震事件分布相叠合，以分析两者之间的联系.图 7b、c和d印证了地震事件往往发生在地热梯度带上的观点，此外，雅安—德阳一带地震事件尤为密集，而同为地热梯度带的甘德县—玛沁县地区地震事件尤为稀少，本文针对这一现象，结合其他地球物理探测结果(Huang and Zhao, 2006；覃庆炎，2011；朱迎堂等，2008；陈飞，2017)总结如下.雅安—德阳一带除位于热流值高-低-高梯度转换带之外，还位于龙门山碰撞造山带之上，是松潘块体和四川盆地两构造分界线，该地区地下两侧物性差异明显、界线清晰且延伸较深：在10 km及以上地层，雅安—德阳一带东侧的四川盆地地区为低速-低阻特征，西侧的松潘构造带为高速-高阻特征；在深部约15~50 km之间，松潘块体的低速-低阻特征对应四川盆地的高速-高阻特征.甘德县—玛沁县一带整体位于松潘块体之上，两侧在横向上不存在物性差异.为了进一步总结规律，本文调研了同为地热梯度转换带的临沧—楚雄一带的地球物理资料(Huang et al., 2006；秦维秉，2017；李冉等，2014；胥颐等，2013)，其物性特点为：该区域下方为低速-高阻体，两侧为高速-低阻体，随着深度的加深，地层的物性差异逐渐减弱.

综合以上讨论可以得出：地震事件往往发生在地热梯度带地区，反之地热梯度带地区不一定易发生地震事件，只有当两侧地层在一定深度范围内存在明显物性差异时，地热梯度带才存在密集的地震事件.雅安—德阳一带正是由于其地热梯度带这一特殊位置以及两侧地层明显的物性差异，使得在受到来自深部的热能上升影响与西北部青藏高原的物质流的影响时，在物性过渡带地区应力传递受阻而产生积累，从而引发地震.

5 结论

经过以上的数据分析与讨论，本文得出以下结论：

(1) 本文计算出青藏高原东南缘地壳上下层热导率、地表和地壳内部的热流分布.热导率低值区(藏南—滇西、乡城、康定和德阳等地)与该地区的地层高温对应较好.地表热流整体上呈现西南高、东北低的变化趋势，研究地区热流异常的高低分布与前人研究成果吻合较好，仅有小部分地区(如康定北段显示非热流高值)与前人结论(热流高值)不符.地壳内部热流值随深度的增加而整体降低.

(2) 通过将地壳不同深度处的热流值与地震事件叠加，本文总结出规律：虽然地震事件往往发生在地热梯级带地区，但是，在地热梯度带地区却不一定存在高频的地震事件，只有当附近地层在一定深度范围内存在物性分布差异时才易出现高频的地震事件.因此，地热梯度带为我们识别孕震区域提供了有利视角，结合其他地球物理探测技术，有望成为圈定易震区域的有力约束.

致谢  衷心感谢两位匿名评审专家提出的宝贵意见.