松辽盆地北部地热场深部控制因素研究：基于大地电磁探测的结果

牛璞, 韩江涛，2*， 曾昭发，2， 侯贺晟， 刘立家，2，马国庆，2， 管彦武，2

1 吉林大学地球探测科学与技术学院， 长春 130026 2 自然资源部应用地球物理重点实验室， 长春 130026 3 中国地质调查局中国地质科学院地球深部探测中心， 北京 100037

0 引言

地热能是一种清洁可再生能源，在替代碳氢化合物资源及发电或直接供热方面有重要作用(Kana et al．，2015；Olasolo et al．，2016；陈昌昕等，2019)．自中生代晚期以来松辽盆地由于岩石圈减薄、广泛的火山活动及花岗岩侵入的共同作用(Wang and Li，2018)，地热异常明显增高，盆地赋存有大量地热资源(刘晨璞等，2016；周立岱，2005)，盆地北部更是处于异常高值区．松辽盆地地温场总体特征为：盆地中央无论是地温梯度、地层温度还是大地热流均呈现高值，并大致以环状向外依次降低(图1)(李野，2017)．因此，研究松辽盆地地热场对于今后地热能的利用具有重要意义．

图1 研究区地热背景及大地电磁测点分布(改编自符伟等，2019；孙成城，2019；李野，2017)(a) 松辽盆地大地热流分布图； (b) 大地电磁测点分布．F1：嫩江断裂；F2：大安—德都断裂；F3：孙吴—双辽断裂；F4：哈尔滨—四平断裂；F5：讷河—绥化断裂；F6：滨州断裂； F7：扎赉特—吉林断裂．Fig.1 Geothermal background and MT sits in the study area (modified from Fu et al.，2019；Sun，2019；Li，2017)(a) Geothermal flow map of Songliao Basin; (b) Distribution of MT sites. F1：Nenjiang fault；F2：Da′an-Dedu fault；F3：Sunwu-Shuangliao fault；F4：Harbin-Siping fault；F5：Nehe-Suihua fault；F6：Binzhou fault；F7：Zalite-Jilin fault.

地热系统通常包含“源、通、盖、储”四要素(杨林和姜国庆，2020)．关于松辽盆地地热场的控制因素可以归纳为以下方面：(1)中国东北地区的深部地球物理结果表明，西太平洋板块深俯冲脱水作用引发软流圈上涌(田有等，2019；Zhang et al.，2020)，松辽盆地发生普遍的岩石圈减薄(韩江涛等，2018)，岩石圈拆沉触发的上地幔热上升是盆地形成热异常的主要原因(Kang et al．，2016；Wang and Li，2018)，其中上涌的部分熔融热物质作为盆地地热系统的热源(刘晨璞等，2016；伍小雄，2014；朱焕来，2011)；(2)区域重磁勘探结果揭示出松辽盆地存在大量的基底断裂(孙成城，2019；邢大全，2015；裴明波，2008；张丛，2017)，可以作为沟通深部热源与浅部热储的主要流体通道；(3)盆地地热储集作用与其稳定盖层关系密切，松辽盆地演化过程中，青一段和嫩一、二段两次大范围湖侵期形成了全盆地分布的泥质岩，致密、低渗的厚层泥岩凭借热导率低的特点作为地热系统盖层，为热物质的保存提供有利条件(刘晨璞等，2016)；(4)盆地北部浅层发育有上下2套热储体系，下部热储体系为泉头组三段、四段地层，上部热储体系为青山口组的青二三段和姚家组地层(胡霞和吕建才，2016)．近年来的地质、地球物理工作较好地阐明了松辽盆地物质组成和结构分布，但是地质结构对地热系统的控制作用仍然不清楚．

1 地质背景

松辽盆地处于西伯利亚板块、太平洋板块、华北板块的构造叠合部位(葛荣峰等，2010)，受到多种构造体制作用，具有复杂的地质演化背景．由于古亚洲洋的闭合，松辽地块在石炭纪末与众多微陆块拼合形成东北地区．450～350 Ma 松辽地块进入构造演化相对稳定时期，发育了相对完整的上古生界(余和中，2001；王成文等，2009)．之后盆地经历了晚侏罗世初始张裂阶段、沙河子组沉积断陷阶段、营城组沉积萎缩期、营城期末坳陷阶段、嫩江组沉积构造反转阶段等复杂的构造演化(王燮培和严俊君，1996)，最终形成叠置于古生代基底上的大型中、新生代裂谷型盆地．

复杂的演化历史为盆地形成地温场提供有利条件．盆地演化初期，受太平洋板块俯冲影响，区域拉张应力导致深部地幔隆起，热物质上涌，盆地岩石圈减薄，造成许多地壳断裂(胡望水等，2005)．演化的中后期，随着沉积序列发展，白垩纪沉积地层逐层覆盖于古生代结晶基底上，形成广泛分布于盆地内的泥质岩盖层．贯穿于整个沉积过程的中新生代火山活动，形成的大量岩浆岩，也为地温场的形成提供有利条件．松辽盆地经历了加里东期、海西期、印支期、燕山期、喜山期等多期岩浆活动(张兴洲等，2012)．在松辽盆地北部地区，大范围缺失三叠系，仅有可能出现华力西期和燕山期火山岩和岩浆活动(朱焕来，2011)．盆地基底岩性主要由弱变形和变质显生宙花岗岩、片麻岩及古生代沉积地层(包括砂岩，灰岩，千枚岩，板岩及变质岩石)组成(Wu et al．，2001；Zhou and Wilde，2013；Wang et al．，2014)．片麻岩可以作为矿化水和热水的含水层(Ozen et al．，2010)．花岗岩在所有基底岩性中放射性元素含量是最高的，其产热率也是所有岩石中最大的(朱焕来，2011)．有学者认为片麻岩的存在，证实盆地基底是前寒武纪的．盆地北部徐家围子断陷的火山岩样品中发现古老锆石，证实松辽盆地基底是由前寒武纪结晶岩系与多期古生代—晚中生代早期岩浆岩构成的复合型基底(章凤奇等，2008)．

2 研究方法

2.1 MT数据

本文基于松辽盆地北部覆盖了地热异常变化强烈区域的大地电磁剖面(图1a)，研究了松辽盆地北部深部地热系统的控制要素．剖面长246 km，共采集71个测点，点距约为3.6 km．剖面自大庆市开始，向北经过林甸县、终到依安县附近(图1b)．地质单元上自南向北依次经过古龙断陷、常家围子断陷、林甸断陷和依安断陷．野外数据采集使用加拿大Phoenix公司的V5-2000型大地电磁仪，采用张量方式布极，每个测点测量3个相互正交的磁场分量(Hx，Hy，Hz)和2个相互正交的水平电场分量(Ex，Ey)，平均采集时间约为20 h．对数据的原始时间序列进行快速傅里叶变换、Robust估计(Egbert，1997)等处理手段，获得阻抗张量信息．最终得到周期范围为0.006～1000 s的有效MT数据，剖面部分测点视电阻率相位曲线(图2)显示数据质量较好．

图2 大地电磁剖面部分测点视电阻率相位曲线Fig.2 Apparent resistivity and phase curves of some MT sites along the profile

2.2 维性分析和电性主轴判别

利用相位张量二维偏离度β进行区域维性分析(Booker，2014；Becken and Burkhardt，2004)．相位张量不受浅部异常体的影响(Caldwell et al．，2004)．使用最大相位(Φmax)作为长轴、最小相位(Φmin)作为短轴绘制相位张量椭圆．主轴方向(Φmax)表示导电结构的走向．通常情况下，使用偏离度β来推断地下介质的维性，将大于3°的值解释为指示三维结构，将小于3°的解释为近似二维或准二维结构(Caldwell et al.，2004)．图3将计算出的相位张量椭圆表示为频率-距离的横截面．随着频率的减小，对深部结构更为敏感．椭圆颜色以偏斜值β着色，研究区电性结构整体成层性较好．大于0.01 Hz时，具有较小偏斜值，指示剖面具有较好的二维性．小于0.01 Hz时，呈现三维特性，所以最终使用0.01～100 Hz的数据进行反演．

图3 相位张量椭圆黑色倒三角-MT测点．Fig.3 Phase tensor map Black inverted triangle: MT site.

相位张量分解、GB分解是判别电性主轴的有效方法(蔡军涛等，2010)，全剖面测点0.001～0.01 s、0.01～0.1 s、0.1～1 s、1～10 s、10～100 s、100～1000 s六个频段的电性主轴玫瑰图如图4和图5所示．可以看出，相位张量分解的结果(图4)，除高频段(0.001～0.01 s、0.01～0.1 s)指示电性主轴方向为0°和90°左右外，其余四个频段指示电性主轴方向大致相同．GB分解的结果(图5)与相位张量分解结果类似．将玫瑰图指示方向、相位张量椭圆长轴指示方向，与区域实际地质走向三者相结合，最终确定区域构造走向为NE35°．

图4 各频段构造走向分析结果玫瑰图(相位张量分解)Fig.4 Rose diagrams showing strike analysis results for different frequency bands (phase tensor decomposition)

图5 各频段构造走向分析结果玫瑰图(GB分解)Fig.5 Rose diagrams showing strike analysis results for different frequency bands (GB decomposition)

2.3 二维反演

基于MT-Pioneer平台集成的二维非线性共轭梯度反演算法(Rodi and Mackie，2001)进行二维反演．在考虑地形变化的前提下，进行TE+TM模式反演．维性分析结果(图3)显示，测点在较低频时具有一定三维性，数值模拟计算结果显示三维畸变对TE模式视电阻率数据易产生影响(蔡军涛和陈小斌，2010)．所以反演时，增大TE模式视电阻率、相位本底误差，降低TE模式对整体反演结果的影响(梁宏达等，2017)．

反演参数：网格个数为47×99，频率范围为0.01～100 Hz．反演过程中设置TM模式视电阻率本底误差为10%，相位本底误差为5%，TE模式视电阻率本底误差为80%，相位本底误差为10%．水平平滑控制参数α为1，垂直平滑度指数因子β为0．

设置初始背景电阻率为10 Ωm、100 Ωm、500 Ωm、1000 Ωm、10000 Ωm的均匀半空间，分别进行TE+TM二维反演(图6)，结果表明，不同的初始模型对反演结果的影响有限，100 Ωm的均匀半空间得到的反演结果相较于其他，结果中冗余构造较少，拟合差最小，模型更加符合实际情况，故最终选择初始模型为100 Ωm的均匀半空间．对不同正则化参数τ(1、3、5、10、15、50、100、300、500、1000)进行测试，以找到模型粗糙度与最终RMS之间的最佳平衡．通过绘制L曲线(图7)，认为拐点处τ值为最佳，故选择τ为15的模型．初始拟合差为13.991，经过260次迭代后，最终电阻率模型拟合差为2.219，拟合差随迭代次数变化情况如图8．

图6 不同背景电阻率二维反演结果Fig.6 Two-dimensional inversion results of resistivity in different backgrounds

图7 不同正则化因子对应模型粗糙程度-拟合差Fig.7 Model roughness and RMS values for different regularization factors

图8 RMS随迭代次数变化曲线Fig.8 RMS variation with iteration number

2.4 模型验证

图9为剖面实测数据和TE+TM二维反演模型的视电阻率与相位的拟断面图，虽然反演中给TE模式视电阻率、相位设置了较大本底误差，使得TE模式数据拟合情况较TM模式拟合情况差，但是整体来看，实测数据和反演模型拟合较好，初步证实了二维反演模型的可靠性．TE和TM模式实测数据和响应数据都显示在浅部整体为低电阻率层，大约0.32 Hz以下电阻率明显增高，与反演电阻率模型吻合．

图9 TETM视电阻率与相位的(a，c，e，g)原始数据(b，d，f，h)响应数据拟断面图Fig.9 Pseudosection sections of observed (a，c，e，g) and modeled (b，d，f，h) TETM apparent resistivity and phase

图10 灵敏度测试结果(a) 高导体C1验证后不同测点RMSUpdate； (b) 高导体C2验证后不同测点RMSUpdate； (c) 剖面二维反演模型高导异常灵敏度测试．Fig.10 Sensitivity test results(a) RMSUpdate of different MT sites after high-conductor C1 verification; (b) RMSUpdate of different MT sites after high-conductor C2 verification; (c) High-conductivity anomaly sensitivity test of profile 2D inversion model.

根据本文2D反演结果设计相应理论模型(图11a)，紫色区域代表反演结果中高阻层(10000 Ωm)，橘色区域代表反演结果中低阻部分(10 Ωm)．对理论模型进行正演，得到该模型响应数据，利用这一响应数据进行TE+TM二维反演，反演参数与实测数据反演时相同，得到反演结果如图11b．与实测数据的反演结果(图11c)进行对比，认为异常体结构位置相似，证明本文反演结果的可靠性．

图11 理论模型验证结果(a) 理论模型示意图； (b) 理论模型反演结果； (c) 实测数据反演结果．Fig.11 Theoretical model verification results(a) Schematic diagram of theoretical model; (b) Inversion results of theoretical model; (c) Inversion results of measured data.

3 电性结构分析

剖面由南向北分别为古龙断陷、常家围子断陷、林甸断陷和依安断陷．电阻率模型总体呈现“纵向分层、横向分块”的特点．以林甸为界，两侧呈现明显不同的电性结构．南侧整体呈现“低阻-高阻-低阻”的三元电性结构，北侧为“低阻”的一元电性结构(图12e)．古龙断陷具有“低阻-高阻-低阻”的三元电性结构特征．从地表向下延伸到-5 km为低阻层，电阻率约为1～10 Ωm．低阻层之下存在一向北延伸的地壳尺度的高阻体R1，厚度近10 km，电阻率大于1000 Ωm．高阻体下存在高导体C1，C1电阻率小于10 Ωm．常家围子断陷电性结构与古龙断陷相似，也为“低阻-高阻-低阻”的三元结构．不同之处在于，浅部低阻层厚度随着向北延伸稍有增厚，高阻体R1厚度向北逐渐变薄，高导体C1呈现小幅度隆升．林甸断陷浅部低阻层厚度向北继续增厚，深部电性结构则以中部林甸为界，南北两侧明显不同．南侧电性结构是与古龙断陷、常家围子断陷相似的“三元”电性结构，而北侧不再有地壳尺度的高阻体，且上地幔存在另一高导体C2，埋深约30 km．依安断陷具有“低阻”的一元电性结构．高导体C2向上延伸至浅部低阻层，两者共同构成了20 km厚低阻体．上地幔不再有高导异常存在，整体呈现低电阻率．

图12 研究区剖面二维反演得到电性结果(a，b)东北地区不同深度S波速度异常(来自Guo et al．，2018)； (c) 地形和构造单元；(d)沿剖面大地热流(红点)； (e) 剖面二维反演电阻率模型(T5、Moho面根据孙成城，2019)．Fig.12 Profile electrical results obtained by 2D inversion in the study area(a，b) S-wave velocity anomalies at different depths in Northeast China (from Guo et al．，2018); (c) Topography and tectonic units; (d) Geothermal flow along the profile (red dots); (e) 2D inversion resistivity model (T5 and Moho from Sun，2019).

4 讨论

4.1 松辽盆地北部地区地热场控制因素

本文将从地热控制要素的四个方面对松辽盆地北部地热系统进行分析．

普遍认为松辽盆地主控热源为幔源热(朱焕来，2011；刘晨璞等，2016；Song et al．，2018)，总地表热量中约65%是盆地下方地幔贡献的(Wang and Li，2018)．本文二维电阻率模型揭示盆地北部下地壳、上地幔存在高导体C1、C2(图12)．通常情况下地球深部的高导体有金属硫化物矿体、结晶石墨膜(Li et al．，2020，2003)、含水流体和部分熔融(Wannamaker et al．，2009)等几种成因解释．松辽盆地具有较高的地表热流值(86～42 mW·m-2)(李野，2017)，而硫化物和互连的结晶石墨膜不能长期在高温下保存(Li et al．，2020)，且金属硫化物矿体成矿规模较小，结晶石墨膜通常为薄层，两者都与高导体C1、C2的赋存形态不同，所以排除这两种可能．东北地区三维速度模型显示松辽盆地中央25 km深度处存在低速异常(Guo et al．，2018)，与高导体C1所处位置相吻合(图12b)．综上，认为C1、C2与流体有关，可能是含水流体或部分熔融．盆地内大量的钻探结果显示深部含有丰富的CO2气藏，天然气的He同位素、碳同位素等化学分析都表明了CO2气藏为无机成因幔源特征(杨悦，2019；付晓飞等，2010)．而且在松辽盆地观察到的深部地幔异常通常归因于部分熔融层或热流体(韩江涛等，2018；杨悦，2019)．结合研究区深部动力学背景，在早侏罗纪古太平洋板块开始高速向西俯冲(Wang et al．，2019)，俯冲板内的含水矿物、断层发生再活化引发板内脱水作用，以及形成的大地幔楔内的对流循环共同作用造成了部分熔融物质的上涌(田有等，2019)，上涌的地幔物质作为主要热源，在电性上表现为高导体，故本文认为深部高导体C1、C2可能是从地幔上涌而来的部分熔融热物质，充当盆地地热系统热源．

地壳断裂是地幔热物质上升至盆地内部的主要流体通道(朱焕来，2011；刘晨璞等，2016)．松辽盆地内发育丰富的基底断裂和地壳断裂，众多的北北东向和北西向基底断裂派生出一些次级小断层、小裂隙，为下部热流向上涌入沉积层创造条件．研究区内主要发育两条深大断裂，大安—德都断裂F2，呈NNE 向延伸并纵贯松辽盆地，向东倾，倾角较陡(Yu et al．，2015)．与之交汇的滨洲断裂F6，属于地壳断裂，倾向北东，呈北西向延伸横穿松辽盆地(孙成城，2019)，是研究区内重要的导热通道．

地热系统中的盖层对于深部热物质的保存有着重要的作用．松辽盆地具有稳定沉积的白垩纪地层，盆地未经历明显的沉积和变形(Wang et al．，2016)．地震反射结果显示盆地总沉积厚度达到8 km(Li and Liu，2015)，电阻率模型浅部低电阻率层的分布特征与地震结果吻合(图12e)，可能代表了盆地内沉积地层．由于白垩纪时期两次大规模湖侵作用，导致盆地内沉积地层主要表现为湖泊-三角洲体系为特征的泥岩沉积层(刘晨璞等，2016)，具有较厚、致密、低渗特点，可作为盆地地热系统盖层，为深部热物质提供良好的保存条件．

热储指的是埋藏于地下、具有有效孔隙和渗透性的地层、岩体或构造带，是热流富集、储存的重要空间(朱焕来，2011)．电阻率模型揭示剖面南部存在地壳尺度高阻体R1，R1向北延伸至林甸附近．松辽盆地内深反射地震结果揭示断陷层(T5反射轴)之下为盆地基底(符伟等，2019)(图12e)，盆地基底主要包括前寒武纪结晶岩体(片麻岩、花岗岩等)和古生界沉积地层(Wang et al．，2014；Zhang et al．，2017；章凤奇等，2008)，对比Lichoro等(2019)在肯尼亚北部火山裂谷带的MT测量结果发现前寒武结晶岩体具有高电阻率特征，因此推测R1为前寒武纪结晶岩体的反映．该结晶岩体内片麻岩具有储热的能力(Ozen et al．，2010)，可以为热能的储存提供有利的空间．

4.2 地质结构与地表热流特征联系

在东北地区，认为松辽盆地岩石圈的扩张和变薄程度是最大的(Meng，2003)，莫霍面(Moho)深度通常在26～32 km之间(Zhang et al．，2020)．此外，中国东北部的岩石圈-软流圈边界(LAB)表现出与Moho面相似的结构特征，松辽盆地LAB埋深最浅约70～100 km(Guo et al．，2014)．在相同的时间，地壳减薄处向上传导热流值较地壳厚地方大．在本文研究区范围内深反射地震结果揭示Moho面埋深变化不大(图12e)，所以在研究区内进一步排除了Moho面变化对于地热分布的影响．从全球热流数据库 (http:∥www.heatflow.und.edu/data.html)以及已发表的中国大陆地区大地热流数据(汪集旸和黄少鹏，1990；胡圣标等，2001；姜光政等，2016)收集到本文剖面所穿过地区的大地热流数据，并通过插值绘制曲线(图12d)．可以看出大地热流的变化与电性结构的变化较好地吻合．大地热流从常家围子断陷开始向北逐渐降低，深部高阻结晶岩体R1厚度也相应地逐渐减薄，直至林甸北侧大地热流最低，此处深部不存在结晶岩体．所以本文认为松辽盆地北部地热场变化主要受地下不同结构的控制．林甸南侧深部存在较厚前寒武纪结晶岩体，为深部热物质的储存提供空间，而林甸北侧深部缺少储存热量的结晶岩体，大地热流明显降低．

4.3 松辽盆地北部地热系统模型

根据本文揭示的电阻率模型，结合深反射地震资料提出松辽盆地北部地热系统模型．松辽盆地北部属于水热型地热资源，主要是由大安—德都、滨州壳断裂等深大断裂控制地热系统形成及温泉分布．林甸地区位于两大断裂交汇处，可直接通过断裂将热流从深部热源向浅部热储传递，存在有众多的温泉．研究区的主要热源是，在晚白垩世古太平洋板块的深俯冲运动后，俯冲板片脱水作用引发的上涌的部分熔融热物质．随后的岩石圈拆沉作用使得热物质向上传递至地壳．盆地发育的深大断裂及派生的次级断裂作为地热系统的流体通道将热流从地壳进一步传递至浅部沉积层及地表．全盆地覆盖的泥质岩是盆地地热系统的盖层．林甸南部的结晶基底具有储热导热作用，将热物质从深部热源传递至上覆沉积层，并储存在浅部热储有利区(胡霞和吕建才，2016)(图13)．

图13 松辽盆地北部地热系统模型示意图Fig.13 Geothermal model of the northern Songliao Basin

5 结论

电阻率模型揭示了以林甸为界的两个不同电性结构．林甸以南的地区地壳存在大规模前寒武纪结晶基底．结晶基底具有储热和导热的作用，使得该地区热储含量和大地热流都为高值．林甸以北地区由于缺少储集热量的结晶基底，热异常逐渐降低．水热型资源富集的林甸地区，位于结晶基底末端，同时还位于两大断裂交汇处，可直接通过交汇的两大基底断裂将热量直接传递到浅地表．

致谢感谢审稿专家提出的修改意见．