遥感影像地表温度反演与地热资源预测——以石家庄地区为例

辛 磊 刘新星 张 斌

1. 中国地质科学院地质力学研究所, 北京 100081;2. 自然资源部古地磁与古构造重建重点实验室, 北京 100081;3. 中国地质大学 (北京), 北京 100083;4. 河北省战略性关键矿产资源重点实验室, 河北 石家庄 050031;5. 河北地质大学地球科学学院, 河北 石家庄 050031;6. 河南省航空物探遥感中心, 河南 郑州 450053

0 引言

地热资源是一种可再生清洁能源, 将在中国能源结构调整中发挥重要作用 (王贵玲等,2020)。 常规地热资源勘查以调查近地表地质现象为主, 辅助以钻探、 物探等方法, 反演、 推断地球深部地质条件, 了解地热资源赋存状况。 然而部分地区受第四系掩盖影响, 利用物探、 钻探等方法来勘查、 评价地热资源, 具有投资大、 风险高等问题 (张中言, 2010)。 因此, 需要借助更加经济、 实用的技术手段, 与其他技术方法相配合,为进一步地热资源勘查提供靶区。

热红外遥感作为一种空间信息探测技术, 其波段数据具有信息量大、 检测精度高、 直观、 形象、 受通行条件限制小等优点, 被广泛应用到多个学科领域研究中 (周彦儒, 1998), 尤其在地热资源勘查领域已取得显著的社会经济效益 (张佩民和张金良, 2006; 连胤卓, 2007; 张中言,2010; 熊永柱等, 2016; 闫佰忠等, 2017)。 在热红外遥感技术中, 相较于单通道算法 (Jimenez-Munoz et al. , 2009) 和基于影像的反演算法 (丁凤和徐涵秋, 2008), 单窗算法 (覃志豪等,2001) 具有较强的实用性 (黄妙芬等, 2006;Wang et al. , 2019), 因此文章选用单窗算法开展研究。

华北地区地热资源的主要类型为沉积盆地型,盆地内热储层多、 厚度大且分布较广, 赋存有大量中低温热水资源 (王贵玲等, 2020)。 而石家庄地区断裂构造发育、 水资源丰富, 地热资源有待进一步查明、 开发利用。 此次研究在河北省石家庄市及周边地区, 以Landsat 8 卫星影像和相关气象数据为数据源, 利用单窗算法反演研究区地表温度, 引入ASTER 夜间热红外影像, 将反演的地表温度图与夜间热红外影像对比 (Eneva and Coolbaugh, 2009; Mia et al. , 2018; Wang et al. ,2019), 得出昼夜异常区, 再将异常区与地质资料复合验证, 研究地热异常的来源及控制因素, 进而圈定地热资源预测靶区。

1 研究区概况

1.1 地理概况

研究区位于冀中南地区, 西部依托太行山,东部连接河北平原, 滹沱河从研究区北部穿过。其地理坐标范围: 东经113°45′～115°19′、 北纬37°30′～38°26′, 面积约为12.9×103km2(图1)。研究区地势总体呈西高东低阶梯状展布, 主要有低山、 丘陵和山前平原三种地貌类型, 属季风性暖温带半湿润大陆性气候, 春季多风沙, 夏季秋季高温多雨, 冬季寒冷干燥。

1.2 地质概况

研究区位于华北板块 (中朝准地台) 中南部(图1c), 太行山山前深大断裂以西属燕山台褶带的东南边 (Ⅱ) 和山西断隆 (Ⅱ) 的东部边缘,以东属华北断拗 (Ⅱ) 的西部边缘 (胡君春和郭纯青, 2007; 张亚东等, 2011)。

太行山山前断裂带位于太行山山脉和华北平原过渡带, 是太行山隆起与华北凹陷的分界断裂,其纵贯全区, 对研究区内沉积作用和构造演化具有明显的控制作用。

太行山山前断裂以西为相对上升盘, 称为西部山区, 新生代地层和沉积物覆盖于太古宙至中生代地质体之上。 新太古代变质岩出露在研究区西北部, 崇礼岩群岩性组合为透辉斜长角闪岩、黑云角闪斜长片麻岩、 角闪变粒岩; 阜平岩群元坊岩组岩性组合为黑云斜长片麻岩、 黑云斜长变粒岩、 斜长角闪岩。 元古代地层集中出露在井陉县东南部, 下元古界为变质结晶基底, 岩性为角闪石片岩、 绿帘石岩; 中元古界为盖层, 包括串岭沟组页岩、 大红峪组石英砂岩和高于庄组碳酸盐岩。 下古生界出露在井陉矿区一带, 寒武系包括馒头组紫红色页岩、 炒米店组粒屑灰岩、 崮山组薄层灰岩和三山子组白云岩; 奥陶系包括北庵庄组泥晶灰岩、 峰峰组灰岩。 区域上缺少晚古生代和中生代地层。

太行山山前断裂以东为相对强烈下降盘, 也称为华北平原, 第四纪堆积物分布广且厚度较大,多在300 ～600 m 之间。 第四系下伏地层有新近系砂岩、 砂砾岩和奥陶系、 寒武系碳酸盐岩。

中酸性侵入岩均发生高级变质作用, 太古代奥长花岗质片麻岩出露在岗南水库西部, 英云闪长质片麻岩、 正长花岗质片麻岩分布在赞皇北部,古元古代变质正长花岗岩在许亭一带, 以许亭岩体为代表。

1.3 地热模式

地热田主要由盖层、 热源、 热储及热水通道组成。 华北平原盖层为地表浅部比较松散、 沉积物仍未固结成岩、 热导率较低且能够阻止热量散失的岩层。 热源为地壳内侵入岩供热或花岗岩中放射性元素衰变产热。 热储为具有一定透水、 储水能力的岩层。 构造断裂带可作为一种导水导热通道, 流体深循环热对流现象常沿断裂带出现,并且断裂交汇部位可作为点式温泉通道 (Tannock等, 2019), 这些地质认识是地热遥感预测的依据。

在华北地区依据地幔柱活动影响和区域构造控制, 地热可分为华北断陷区深部地热和华北内陆造山带热泉地热 (牛树银等, 2001)。 燕山运动以来, 华北平原受地幔亚热柱强烈上涌影响, 地壳形成一系列相间排列的二级地堑、 地垒构造和三级凹凸构造, 由于多期次岩浆、 岩脉侵入作用,在岩浆上部或局部隆起区、 堑垒边界密集断裂带处易形成地热田; 其外围受幔枝构造作用和岩浆侵入作用影响, 造山带内发育密集断裂, 断裂区可以形成地热泉循环系统 (王钧等, 1983; 牛树银等, 2001)。 基于不同构造成因和热传递方式,地热成藏模式可分两类, 京津冀隆起山地对流型地热资源成藏地热模式和京津冀沉积盆地传导型地热资源成藏地热模式 (图2), 模型详细介绍见王贵玲等 (2017), 更多地热能传输机制可参考罗文行等 (2020)。

图2 京津冀地区地热资源成藏模式 (据王贵玲等, 2017 修改)Fig.2 Accumulation model of the geothermal resources in the Beijing-Tianjin-Hebei region (modified after Wang et al. , 2017)

地表温度是热红外遥感地热预测不可或缺的参数, 可以量化异常与非异常之间的数量关系。一些学者研究 (许军强等, 2008; 张中言, 2010;闫佰忠等, 2017) 表明, 地热田是地壳中热量在断裂破碎带中运移、 循环而形成的。 值得指出的是, 正是由于地壳深部热量部分沿断裂构造传递到地表, 即热异常顺构造行迹产出, 才可以被遥感技术探测、 识别。

因此, 需要从两方面展开研究: 基于遥感影像的热红外波段反演温度异常区域和利用遥感影像解译研究区内断裂构造圈定地热资源远景区。

2 基于热红外波段地表温度反演

2.1 Landsat 8 和ASTER 数据特点

Landsat 8 数据: OLI (陆地成像仪) 有9 个波段, 成像宽幅为185 km×185 km。 TIRS (热红外传感器) 有2 个波段, 空间分辨率为100 m (USGS,2013)。 由于Landsat 8 TIRS 的Band 11 辐射校正偏差较大 (USGS, 2014, 2016), 因此选用Band 10进行研究。 波段基本情况见表1、 表2。

选用了2 幅可完全覆盖研究区的Landsat 8 卫星影像数据 (影像基本信息见表3、 表4) 进行镶嵌、 裁剪及其他预处理操作。

2 幅ASTER 卫星影像数据的Band 13 影像, 经预处理后, 用以降低日间地表异常温度解释的多解性。 ASTER Band 13 影像波段范围 10.25 ～10.95 μm; 量化等级12 bits; 空间分辨率90 m;ASTER 数据获取时间分别为2013 年1 月3 日14 ∶21 ∶35.186 (GTM) 和14 ∶21 ∶44.026 (GTM)。

表1 TIRS (热红外传感器)Table 1 TIRS (Thermal Infrared Sensor)

表2 OLI (陆地成像仪)Table 2 OLI (Operational Land Imager)

2.2 地表温度反演

2.2.1 地表亮度温度反演公式

覃志豪等 (2001) 单窗算法的优点在于: 地表温度反演结果好 (黄妙芬等, 2006; Wang et al. , 2019), 且仅需要3 个基本参数, 分别为地表比辐射率、 大气透射率和大气平均作用温度。 单窗算法公式如下:

其中,Ts是地表温度, 适用范围为273～343 K;a、b为回归系数, 分别为a=-67.35535,b=0.458608;ε是地表比辐射率;τ是大气透射率;Tb是热红外波段的像元亮度温度, K;C和D为参数;Ta为大气平均作用温度。

热辐射强度与Landsat 8 影像DN 值之间的关系可表示为公式 (4):

其中,Lb为传感器接收到的辐射强度 (W·m-2·Sr-1·μm-1),Qdn为Band 10 影像亮度值。 通过辐射亮度公式 (5) 计算亮度温度, 其公式为:

表3 影像1Table 3 Image 1

表4 影像2Table 4 Image 2

公式中,K1和K2为发射前预设的常量。 Landsat 8 TIRS 影像元文件给出, Band 10 的K1=774.89 W·m-2·Sr-1·μm-1,K2= 1321.08 W·m-2·Sr-1·μm-1(覃志豪等, 2001; 张中言, 2010)。

2.2.2 确定大气透射率

大气透射率受多种气象因素共同影响, 例如大气压力、 温度、 气溶胶含量等, 所以计算准确的大气透射率参数不仅复杂, 而且需要详细的气象数据 (连胤卓, 2007; 张中言, 2010)。

此次研究使用的气象数据来自全球天气精准预报网: https: / /www. wunderground. com/history/daily/cn/shijiazhuang/ZBSJ/date/2015-3-6, 数据如表5。

表5 2015 年3 月6 日石家庄正定气象数据Table 5 Meteorological data of Zhengding on March 6, 2015

获取星-地同步观测大气资料后, 大气透射率可以通过辐射传输软件MODTRAN 模拟得到, 并且韦玉春等 (2015) 指出, 在 http: / /atmcorr.gsfc. nasa. gov/页面中输入成像时间及中心经纬度就可估算大气透射率参数。 依据2015 年3 月6 日格林尼治标准时间 (GTM) 3 ∶00 的Landsat 8 卫星影像参数以及同时刻气象数据, 求得大气透射率参数τ =0.96 (表6)。

2.2.3 确定大气平均作用温度

大气平均作用温度主要取决于气温和大气状态。影像数据获取时间为2015 年3 月6 日, 研究区中心纬度约为38°, 因此选择使用覃志豪等 (2003)总结的中纬度冬季平均大气作用温度Ta估算方程:

其中,T0为卫星过境时的地面 (一般为距地面2 m高处) 气温, 单位均为K。

2.2.4 地表比辐射率的确定

丁凤和徐涵秋 (2008) 总结出以下地表比辐射率计算方法: 将遥感影像中地物类型分为水体类、 城镇类和自然表面类3 种类型, 其中水体比辐射率参数为0.995, 自然表面和城镇比辐射率参数需要分别依靠公式 (7) 和 (8) 进行计算求得:

其中,ε地表和ε城镇分别代表自然地表比辐射率和城镇比辐射率,Pv为植被覆盖度, 计算方法如下:

表6 大气透射率估算模型Table 6 Atmospheric transmittance estimation model

式中,NDVI为归一化差异植被指数, 取NDVIv=0.70 和NDVIs=0.05, 当某个像元的NDVI大于0.70 时,Pv取值为1; 当NDVI小于0.05,Pv取值为0。

2.2.5 地表反演温度结果验证

基于上述原理, 经过波段运算求得研究区地表反演温度结果(图3)。 通过像元统计, 可知其平均温度约为17 ℃, 且温度主要集中在12～23 ℃。

图3 地表反演温度图Fig.3 Land surface temperature retrieval map

在GTM 3 ∶00 时, 西部山区裸露岩石热惯量小, 温度上升幅度大, 其温度显著高于研究区平均温度; 东部平原受季节和气候影响, 滹沱河与磁河流量较少, 河床裸露, 温度上升幅度明显,高于平均温度3 ℃; 由于水体热惯量大, 相同辐射条件下研究区内岗南水库和黄壁庄水库温度相对较低, 约为7 ℃。

影像采集时间属于北方冬小麦返青期初期,小麦发黄、 矮、 分布稀疏、 土壤覆盖度低, 导致城镇周边的农田区温度上升不显著, 仅高于平均温度1 ℃, 高于城市气象观测站点实测温度8 ℃。受居民正常生活影响, 城市、 村镇内空气湿度相对高于郊区、 农田, 因此居民区比农田区温度上升幅度低3 ℃。

根据以上现象可知, 地表温度反演结果符合真实情况。

3 遥感地质构造解译

3.1 构造解译

地热资源分布与构造之间存在密切关系 (牛树银等, 2001; 连胤卓, 2007; 许军强等, 2008;张中言, 2010; 王贵玲等, 2017; 闫佰忠等,2017), 需利用遥感影像解译研究区线环构造。 将研究区的Landsat 8 OLI 多光谱影像分别进行线性拉伸、 反差增强和彩色合成等处理, 并依据地质构造直接解译标志、 间接解译标志进行线性构造解译。

解译结果 (图4) 显示, 线性构造主要分布在太行山山前断裂以西。 东部平原与西部山区呈突变式接触, 可推断山区与平原属于断层接触, 该断裂在石家庄以北主要以北北东向展布, 以南为近南北向展布, 中间呈过渡现象。 西部山区解译构造主要为北东、 北北东、 北西和北西西向线性构造, 少数线性构造为近东西向展布。 在平山—井陉—鹿泉一带解译出4 个环形构造。

图4 构造解译结果 (754 彩色合成)Fig.4 Structural interpretation results and previous research results (color composite of bands 7, 5, and 4)

3.2 物探资料对遥感解译结果印证

东部平原广泛分布的第四系对该地区地热资源预测具有阻碍作用, 使得遥感解译该区域地质演化所形成的线性构造较为困难, 然而剩余重力异常反映地壳浅部岩石密度和厚度的变化, 在剩余重力异常图中重力异常带轴向错动或局部等值线扭曲位置处可以推断断裂构造存在。 孙会玲(2019) 认为大比例尺地球物理数据在识别盖层及断裂构造方面有显著优势。 研究区剩余重力异常以北东向展布为主 (图5), 高低相间排列分布反映了深部基底局部的隆起与凹陷, 并且北东向展布的剩余重力高值异常带存在北西向局部错动现象, 可推测石家庄地区北北东向构造格架被北西向断裂错动。 石家庄西部山区北东、 北北东、 北西, 以及近东西向断裂 (＞10 km 的断裂) 的空间分布规律与剩余重力反映的构造格架基本吻合(图5)。 受剩余重力异常数据的精度限制, 无法对遥感解译结果中区域构造之间的次级构造和局部小断裂进行有效验证。 太行山山前断裂以东的平原地区, 第四纪沉积物覆盖严重, 目前无法有效解译平原区构造, 因此需要结合一些学者的研究成果进行综合分析。 解译结果 (图4) 显示, 研究区内以北东、 北西和北西西向断裂分布为主, 北西和北西西向断裂错断北东向断裂, 在研究区西南山区分布近南北向断裂。 平原地区断裂分布在重力异常带等值线扭曲部位 (即正负异常交替带上) 或异常带错断部位 (图5), 可将其理解为控制地壳浅部构造格架的隐伏断裂。

图5 研究区剩余重力异常图Fig.5 Residual gravity anomaly of the study area

4 地温异常对比讨论及地热预测

4.1 识别、 评价温度异常区

众多学者 (Eneva and Coolbaugh, 2009; Mia et al. , 2018; Wang et al. , 2019) 研究表明, 遥感技术地热预测时引入夜间热红外数据, 可降低白天地表异常温度的多解性。 因此, 文中将地表温度图与2013 年1 月3 日的夜间影像进行对比分析, 寻找昼夜温度异常区。 将异常区与地理图关联, 可将地表温度异常区分为两种: 非人为地温异常和人为地温异常 (市、 区、 县和镇)。 人为地温异常区域受建筑物、 城市效应等因素的影响,对这些区域做地热异常预测需要考虑的因素更为复杂, 因而非人为地温异常区是此次研究重点。将非人为地温异常区与地质、 地球物理资料对比分析, 进而得到存在地热的高潜力区域。

4.2 山区寻找地热

依据地热资源成藏模式 (图2), 太行山山前断裂带以西的山区地热资源的分布和热循环受断裂构造控制。 已知研究区中有1 处已开发利用的地热资源——平山县温塘地热田 (胡君春和郭纯青,2007), 现有开发利用地热井8 眼, 地热井开采深度一般为80 ～100 m, 单井涌水量30 ～60 m3/h, 热水水温45 ～73 ℃。 张博 (2015) 对温塘地热田地质特征的研究表明, 温塘地热田受断裂构造控制,北东向断裂破碎带 (倾向南西) 纵穿地热田, 北西西向断裂截切北东向断裂, 是断裂构造破碎带型热储; 北东断裂作为储水构造, 北西西断裂为导水导热构造, 渗水沿断裂带循环过程中不断吸收热源的热量, 在断裂交汇位置处以热水形式侵位于近地表, 承压热水水位埋深14.0 ～20.7 m; 由于温塘地热田上覆第四系松散沉积物厚度小, 为5 ～20 m, 隔温性能差, 属于中低温半封闭式山地对流型地热系统。 韩亚彬和张崇山 (2020) 通过地球物理数据推断温塘地区存在隐伏断裂和中酸性隐伏岩体。 根据温塘地热田遥感影像特征及相关研究显示, 白天地表温度具有相对较高的值,即反演温度平均约为17 ℃, 高于背景温度约3 ℃。该区域夜间影像相对于周围地形、 地物具有较高像元亮度值 (图6)。

根据平山温塘地热田影像特征和地质特征对比结果, 发现研究区内存在1 处与其异常特征相似的地温异常区域——平山县寺沟村区域 (图7)。该区白天平均地表反演温度约为20 ℃, 高于背景值3 ℃, 夜间热红外影像像元亮度值相对高于邻近区域, 北东向解译线性断裂切穿寺沟村区域, 剩余重力异常图中正负异常梯度带可印证该断裂存在。

寺沟村地区断裂构造发育, 可作为良好的导水导热通道控制地壳深部热量沿断裂构造循环至地表, 进而体现在遥感影像上。 该区域汇水面积大且距离岗南水库近, 有利于水源入渗补给, 渗水通过断裂垂直入渗参与深循环, 补给地热田。

图6 平山县温塘镇温泉异常特征Fig.6 Anomaly characteristics of the hot springs in Wentang Town, Pingshan County

图7 平山县寺沟村异常特征Fig.7 Anomaly characteristics of Sigou Village, Pingshan County

4.3 平原地区寻找地热

华北平原内地热流体通过深大断裂使地表水循环至地壳深部, 加热后沿断裂运移至基岩隆起区 (王钧等, 1983)。 地热钻孔数据 (胡君春和郭纯青, 2007) 表明, 石家庄市平原地热探测井深98.2 ～5435 m 不等, 但主要分布在700 ～3000 m 深度范围。 平原地区新生界沉积物覆盖较厚, 有大面积农田覆盖, 以及受其他人为因素影响, 该区域昼、 夜遥感影像地热资源预测可行性不高。

对剩余重力异常图进行构造解译 (图8) 可知: 藁城—无极一带、 辛集以西的马于—换马店一带均为剩余重力正异常, 被周围剩余重力负异常包围, 对应深部隆起构造区域, 这些区域符合京津冀沉积盆地传导性地热成藏模式。 藁城—无极隆起和马于—换马店隆起区位于华北沉积盆地边部, 第四系松散沉积物厚度大, 相对于密度较大的古生界、 中生界和新生界砂岩、 含砾砂岩以及灰岩等岩层, 其热导率低, 可作为隔热层。 隔热层使来自地下深部的热量不易外散而聚集在盖层底部的高孔隙度、 高热导率岩石地层之中。 合理采用热储回灌系统和工艺流程 (阮传侠等,2017) 开发地热资源, 并持续监测、 应对可能因地热资源开采所引发的地质灾害 (甘浩男等,2020), 可以有效保证地热资源可持续开发利用。更多关于京津冀地球物理数据及解译结果见张亚东等 (2011) 和方菲 (2020)。

京津冀地热成藏模式在地热资源预测过程中起着理论支撑作用, 以断裂为通道或纽带, 将地壳浅部的隆起和凹陷构造分别产生的热源场与地表温度相关联, 相互论证, 得出地热远景预测结果 (图8)。 在太行山山前断裂以西的平山温塘地热田处, 总结勘查标志, 并基于勘查标志进行类比, 圈出1 个地热远景区, 即平山县寺沟村一带;在东部平原地区, 结合京津冀地热成藏模式圈定藁城—无极一带、 马于—换马店一带共2 处地热预测远景区, 可为地热资源勘查、 勘探、 定量化评价提供依据。

图8 研究区构造解译与远景预测图 (据张亚东等, 2011; 方菲, 2020 修改)Fig.8 Structural interpretation results and geothermal resource potential of the study area (modified after Zhang et al., 2011; Fang, 2020)

5 结论

基于Landsat 8 数据Band10 影像, 利用单窗算法对石家庄地区进行地表温度反演, 并将地表温度反演结果与夜间热红外数据、 人文地理数据进行对比分析, 排除了人为活动影响, 降低了地表异常温度的多解性。 在此基础上, 综合线性构造解译结果、 剩余重力异常数据及地热成藏模型等,对研究区地热异常进行综合分析, 预测了平山县寺沟村一带、 藁城—无极一带、 马于—换马店一带共3 处地热远景区。

研究表明地表温度反演结果符合实际情况,可以在一定程度上指示研究区地热异常分布。 同时, 在空间上地热异常分布受主要断裂构造控制,可以被剩余重力异常和研究区地热成藏模型印证、解释, 表明地热预测结果具有合理性。

虽然采用热红外遥感技术、 构造解译、 剩余重力异常数据对研究区地热资源进行了深入研究及预测, 但是受天气、 影像数据质量和数据精度等因素影响, 使得预测结果存在准确性、 精确性不足等问题。 因此, 在后续研究中需要利用精度更高的数据, 同时利用更多遥感相关技术对地热资源进行远景预测, 例如: 地热相关的蚀变矿物信息提取、 应用于山区的地形坡向校正等, 实现多元数据综合约束, 以取得更准确、 更精确的预测结果。