Landsat ETM+遥感数据在寻找临沂地热中的应用

刘效才，江泳，刘孝阳，刘琨

(1.山东省第七地质矿产勘察院，山东 临沂　276000；2.临沂市国土资源局，山东 临沂　276000；3.中国地质大学(北京)土地科学技术学院，北京　100083)

Landsat ETM+遥感数据在寻找临沂地热中的应用

刘效才1，江泳2，刘孝阳3，刘琨1

(1.山东省第七地质矿产勘察院，山东 临沂276000；2.临沂市国土资源局，山东 临沂276000；3.中国地质大学(北京)土地科学技术学院，北京100083)

沂沭断裂带内蕴藏着丰富的地热资源，以landsat ETM+为主要遥感数据，利用辐射传输方程反演出整个临沂地区地表温度，对灰色温度图像进行交互式拉伸突显出33℃及以上温度异常区，结合ETM+真彩色图像以及Google影像，排除城镇建筑、农村居民点、河岸滩、裸露岩体等异常区域，将所得异常结果与地热地质背景图进行叠加比对，并从区域地质角度对其进行分析，以期找到地热潜伏区。

Landsat ETM+；温度反演；地热；临沂

引文格式：刘效才，江泳，刘孝阳，等.Landsat ETM+遥感数据在寻找临沂地热中的应用[J].山东国土资源，2015，31(5)：76-81.LIU Xiaocai，LIU Xiaoyang, JIANG Yong，etc. Application of Landsat ETM+ Remote Sensing Data in Searching Geothermal Field in Linyi City[J].Shandong Land and Resources,2015，31(5)：76-81.

0　引言

地热是蕴藏在地球内部的一种自然热能，其热量主要来自地球内部放射性元素的蜕变。从能源的可持续发展战略看，地热资源是一种发展前景非常好的绿色能源，它分布广泛而且性质较稳定，与煤、石油、天然气等传统能源相比，具有廉价、洁净和可持续利用等优点。据估算，地球蕴藏的地热能约14.5×10J，折合4948万亿t标准煤，而全球煤的埋藏量仅为1万亿t[1]。加上近年来全球气候变暖，雾霾天气频发等环境污染问题日益突出，地热资源开发利用的巨大前景备受各国重视。科学开发利用地热能能够实现经济效益社会效益和环境效益的统一，有着各方面的重要意义[2]。

比较成熟的地热资源的探测方法有地质学方法、地球物理勘探方法以及地热异常现象调查法、地球化学法等[3]。但这些传统方法周期长，资金投入大，如果对探测区地热宏观环境条件不明，就会出现投资风险大或盲目开采的情况，对地热资源的可持续发展是一种很大的危害[4]。随着科学技术的不断发展，遥感技术的不断成熟，其在地热资源探测与调查中的优越性越来越明显：可快速、经济、大区域的定量研究热量资源和热场分布变化。郝艳梅等于2001年探索了遥感技术在小区域应用的方法和途径,揭示了遥感技术的广阔应用前景[5]；乔玉良等于2002年利用TM6作为数据源对快速寻找地下热水进行了探究[6]；杨波等于2003年利用遥感技术发现了腾冲西南地区地热富集的规律性，提出了地热田的遥感影像模式[7]；章方等于2012年利用Landsat 7 ETM+遥感影像对长春地区地表温度异常进行了遥感分析[8]。

1　研究区域与数据资料

1.1研究区域概况

临沂市近黄海，东连日照，西接枣庄、济宁、泰安，北靠淄博、潍坊，南邻江苏。南北最大长距228km，东西最大宽度161km，总面积17191.2km2，地势西北高东南低，是山东省面积最大的地级市。临沂市地质构造复杂，地层发育比较齐全，从太古宇至新生界，除上奥陶统、志留系、泥盆系、下石炭统、三叠系及中、下侏罗统、古新统、新近系缺失以外，其他各期地层都有发育。主要构造以郯庐断裂带(境内称沂沭断裂带)为主。沂沭断裂带经郯城北向延伸，纵贯临沂市，以断裂为界，临西为鲁西台背斜的一部分，属华北地台范畴，临东为胶南隆起的一部分，属扬子大陆块范畴。

沂沭断裂带由4条主干断裂组成，自西向东依次为鄌郚-葛沟断裂、沂水-汤头断裂、安丘-莒县断裂及昌邑-大店断裂[9](图1)。该断裂带不但控制了地层展布，岩浆活动也沿断裂侵入，同时把深部热液沿断裂向浅部运移，为地热形成提供了通道和储存空间[10]。复杂的地质构造,为区内地热资源形成提供了地质条件,是地热资源较丰富地区[11]，临沂也因此被国土资源部授予“中国温泉之城”的美称。目前已发现多处温泉，如汤头地热田、铜井地热田、松山地热田等。由于近年来勘查工作一直未能取得突破性进展，加上温泉长期过量开采，多处温泉出现不同程度的枯竭，因此方向正确的温泉勘查工作迫在眉睫。

1—第四系；2—新近系；3—古近系；4—白垩系；5—奥陶系；6—寒武系；7—沂水岩群；8—花岗闪长岩；9—闪长玢岩；10—二长花岗岩；11—韧性剪切带；F1—甘霖断裂；F2—燕甘断裂；F3—蒙山断裂；F4—新泰-垛庄断裂；F5—铜冶店-孙祖断裂；F6—相邸断裂图1　临沂市地热地质背景图巩贵仁等修改)

1.2数据资料

该文所用遥感数据为2003年5月11日的Landsat7 ETM+遥感影像，卫星过境时间为上午10点31分，图像质量较好，清晰无云。所用气象统计资料来源于中国气象局。

2　地表温度反演

2.1辐射定标

ETM+数据是以灰度值来表示的，灰度值并非实际的反射或者辐射强度值，但它们之间存在着一定的换算关系，在进行温度反演之前，必须要将灰度值转化为传感器接受到的辐射强度值，即所谓的辐射定标。对于ETM+数据，辐射定标的转换公式可以简化为

L(λ)=gain×QDN+offset

式中：QDN为像元灰度值；gain和offset分别为图像的增益和偏置，是数据头文件中提供的定标系数。

2.2地表比辐射率计算

物体的比辐射率是物体向外辐射电磁波的能力表征，这一基本参数主要取决于地表物质结构和观测波段，与观测角度等因素也存在关系。ETM+热红外波段分辨率为60m，大都为混合像元，在大尺度上对比辐射率精确测量的难度很大，目前一般采用间接的方法来获取。地球表面不同区域的地表结构虽然很复杂，但从卫星像元的尺度来看,可以大体视作由3种类型构成：水面、城镇和自然表面[12]。根据VandeGriend和Owe的实测数据研究，ε与植被覆盖度存在密切关系。该文将水体像元的比辐射率赋值为0.995，自然表面和城镇像元的比辐射率估算则分别根据下式进行计算：

式中：εsurface和εbuilding分别代表自然表面像元和城镇像元的比辐射率；Fv代表植被覆盖度，指植被冠层的垂直投影面积与土壤总面积之比。采用的是混合像元分解法，将整景影像的地类大致分为水体、植被和建筑，具体的计算公式为：

FV=(NDVI- NDVIS)/(NDVIV-NDVIS)

式中：NDVI为归一化差异植被指数，NDVIv为100%植被覆盖区域的NDVI值，对应研究区域NDVI最大值；NDVIs为100%裸土覆盖区域的NDVI值，对应研究区域NDVI最小值。然而影像不可避免地会受到噪声的影响，可能对NDVI值产生一定的影响，为了保证估算精度，该文利用不同波段彩色合成分别识别出植被与土壤，经过统计分析，得到NDVIv=0.5，NDVIs=0。且有，当某个像元的NDVI大于0.5时，FV取值为1；当NDVI小于0时，FV取值为0。

2.3计算相同温度下黑体的辐射亮度值

在热红外窗口，大气并非是完全透明的。传感器接收的热红外辐射除地表信息以外，还受大气状况的影响。由于热红外波段大气散射一般可以忽略不计，热红外辐射的大气传输方程应考虑大气的吸收和发射[13]。Chandrasekhar[14]把传感器接收的热红外辐射表达为3项之和，即大气上行辐射项，经大气吸收衰减后到达传感器的辐射，大气下行辐射项经地表反射后通过大气传输到传感器的辐射。简化Chandrasekhar的计算表达式得到辐射传输方程：

式中：ε为地表辐射率；TS为地表真实温度；B(TS)为普朗克定律推到得到的黑体在TS的热辐射亮度；τ为大气在热红外波段的透过率。则温度为T的黑体在热红外波段的辐射亮度B(TS)为：

根据遥感图像时间以及中央经纬度等信息，在NASA官网模型中查询式中所需参数值。

2.4反演地表温度

在获取温度为TS的黑体在热红外波段的辐射亮度后，根据普朗克公式的反函数，求得地表真实温度TS=K2/ln(K1/B(TS)+1)。对于ETM+数据，K1=666.09W/(m2·sr·μm)，K2=1282.71K。温度反演结果见图2。

图2　临沂市地表温度图

2.5反演地表温度验证

该文所用辐射传输方程法反演出临沂地区最高温度38.17℃，最低温度16.58℃，大部分地区温度介于20～30℃之间。由于没有卫星过境时的同步气象观测数据，所以不能对地表温度反演的结果进行客观的评价[15]。该文只是采用中国气象局统计的往年5月份气温资料进行简单比对。据统计，临沂地区5月份温度极端最低温度为3℃(1961年)，极端最高温度为38℃(1982年)。从图3可以看出，平均最高温度约为25℃，平均最低温度约为14℃。由于陆地与大气热交换的不平衡性，一般情况下认为白天陆地温度会略高于大气温度，水域温度略低于大气温度，与所反演结果大体吻合。从温度空间分布来看，河流、水库区域温度较低，北部山区和城镇地区温度较高，温度分布符合逻辑，所以认为反演结果准确。

图3　临沂市5月份温度走势图

3　温度异常区分析

3.1温度异常区提取

为了便于提取温度异常信息，该文对灰度值温度图像进行交互式拉伸，突显33℃以上异常区，并用ENVI软件与ETM+真彩色图像进行地理链接。由于ETM影响分辨率为30m，不能很细致地观察地表成分，只能对温度异常区进行粗筛选，排除一些显而易见的人为或地形地貌因素干扰。为了更好地对异常区进行分析，该文将交互式拉伸处理后的图像加载到Global mapper中，对温度异常区进行数字化，圈出共计58个区域，并导出为KMZ矢量格式，利用Google Earth高分辨率和历史图像等优势，对其进一步比较分析。

经认真比对分析，该文将温度异常区划分为两类：A类温度异常区和B类温度异常区。A类温度异常区指由于城镇建筑、村庄、河岸滩、裸露山体等因素影响下呈现的温度异常区，这种温度异常并非地热因素导致。B类温度异常区排除人为或地形地貌等因素对温度的影响，温度异常原因尚不明确，待考证。两者共计58处，其中A类温度异常区39处，B类温度异常区19处，对B类温度异常区进行编号B1-B19。分析结果如图4所示。

图4　温度异常区在Google Earth上的显示

3.2B类温度异常区与临沂市地热资源开发利用规划图的叠加分析

图5为临沂市政府所编制下发《中国温泉之城——山东临沂发展建设规划(2013—2020年)》的地热开发利用规划略图，包含2013年前已经开发的温泉示范区，2013—2015年温泉规划开发区以及2016—2020年温泉规划开发区。其中已开发示范区已完成勘察验证，规划开发区为地热异常区。

图5　临沂市“一轴两翼”地热资源开发利用规划略图

将该文所得B类温度异常区与之叠加分析，共发现6处B类温度异常区完全处于已开发示范区或规划开发区范围内，编号分别为B1，B5，B8，B14，B15，B19，它们显示出的温度异常面积大都较小，呈现为点状分布，应为地热露头或者地热中心区。B4，B10，B12三处温度异常区与温泉已开发示范区或规划开发区有大部分重叠，它们的面积较大，呈现长条带状。其中B10受到北部费县处温泉和区内山体的共同影响，表现出温度异常形状狭长且面积较大；B12处于城市建筑密度较大区域，受其与地热的共同影响，表现出的面积也较大。部分地热区在热红外遥感温度图像上得到一定程度的体现，表现为热异常，但由于受到热源埋深、温度和地表覆盖情况等因素的影响，并非所有地热区都会在地表呈现出热异常，有些地热区未能监测的到。

3.3B类温度异常区与地热地质背景图的叠加分析

排除A类温度异常区和上文中所列与温泉已开发或规划开发区重叠部分的异常区，该文利用图像叠加层功能，将所剩B类温度异常区与临沂市地热地质背景图在Google Earth中进行叠加，从区域地质角度对其进行分析判断地热存在的可能性。经分析，共发现3处B类温度异常区具有存在地热的可能性，且均处于韧性剪切带或者断裂带上，编号分别为B2，B3，B9(图6)。

图6　温度异常区与地质图在Google Earth上的叠加

3.3.1B2号异常区

该区位于蒙阴县东部兴旺庄村—魏家岩路村一带，呈NW—SE向条带状分布，处于铜冶店-孙祖断裂上。

铜冶店-孙祖断裂为NW向旋转构造带的临北帚状构造，全长100km，走向NW，倾向NE；断裂带内构造片岩、糜棱岩、角砾岩发育，断面上斜冲擦痕清晰，断盘上地层挤压牵引明显，不利于地热水流通，也是阻水断裂。N盘(弧外)下降由泰山岩群构成，S盘(弧内)上升由古生界、中生界组成；控制着中、新界凹陷盆地沉积，其特点为北断南超，地层东老西新；由一组多期活动的左行逆断层组成；活动时间为侏罗纪至新近纪，白垩纪活动最强烈。该断裂位于水文地质条件良好，富水性强的沂沭断裂带以西中低山丘陵水文地质区，该区域分布大面积的碳酸盐类裂隙岩溶水，含水层主要由灰岩组成，厚度大，分布广，裂隙岩溶发育，彼此连通，为地下水提供了充裕的赋存空间。铜冶店-孙祖断裂与NNE向断裂交会处有非常好的成热条件，同时由于鄌郚-葛沟断裂为由北向南的疏水断裂，故NNE向断裂与其交会部位为地下水汇集地区，地热水水源丰沛。

3.3.2B3号异常区

该区位于沂南县兰香埠村—司马村一带，呈NNE向条带状分布，处于沂水-汤头断裂带上。

沂水-汤头断裂走向15°，构成苏村-马站地堑的东界和汞丹山凸起的西界，断裂南段出露不好，北段沂水境内出露好。在新生代早期构造处于开放系统，火山活动频繁，在该断裂附近发育一系列派生断裂，在断裂交会处对地热形成极为有利，汤头温泉就在该断裂带上。

3.3.3B9号异常区

该区位于费县盛家岭—安乐窝村一带，呈NE向条带状，处在甘霖断裂处附近；区内地层多以寒武纪海相沉积岩为主，在中生代印支期被中酸性花岗岩沿裂隙或顺层侵入。推断区内中生代侵入岩根部及甘霖断裂提供了地热资源的热源及通道，存在较好的地热勘查前景。

4　结论

该文从热红外遥感和区域地质角度对临沂市地热分布状况进行了简要分析，共有部分温泉已开发区或者规划开发区显示出了温度异常，在遥感温度图像得到了很好的验证。另外发现3处地热潜伏区，对以后临沂地区地热勘测具有一定的导向性作用。在实际应用中，此法与同类技术相比具有以下优势：

(1)利用热红外波段进行地表温度反演，确切的量化出整个临沂地区地表温度，而以往遥感地热相关研究往往是通过波段差值运算和增强处理等步骤定性的区分地表温度高低，对于温度异常区缺乏定量分析。相比之下，通过温度反演这一过程使得研究更加精确。

(2)对于灰色温度图像进行交互式拉伸，凸显出33℃以上区域，而使得其他温度区作为黑色背景值不予显示，使得图像更具有可读性。

(3)该文将33℃以上温度异常区进行了矢量化，并利用Google Earth具有高分辨率和历史图像等优势排除人为因素或地形地貌等引起的温度异常区，把所得到的结果分别与地热规划图和地热地质背景图进行叠加分析，通过2次筛选增加了结果的准确性和可靠性。

但由于缺乏物探、钻探以及磁异常等实际工作，对得出的3处地热潜伏区也只属于初探，尚需等待验证。遥感虽快捷方便，可节省大量人力物力，但毕竟还不能完全代替传统的地热勘测技术。遥感技术必须与其他技术方法相配合,与专业知识及专家思想相结合才能取得比较好的应用效果。

[1]中国国土资源报.国际视野下的地热资源开发利用[J].国土资源，2014，(1):56-57.

[2]吴新雄.科学、清洁、高效、可持续地发展地热能[J].能源，2014，(3):38-39.

[3]曹凯,江南.基于TM6的地热资源的热红外遥感探查模型研究[J].遥感信息,2006，(2):18-21.

[4]乔玉良,吕芳,王栋,等.卫星遥感技术在寻找地下热水中的应用[J].航天返回与遥感,2007,23(4):58-62.

[5]郝艳梅,高新法,张义文.遥感技术在下庄地区地热资源开发中的应用研究[J].国土资源遥感,2001，(1):19-24.

[6]乔玉良.应用红外遥感寻找地下热水初探[J].航天返回与遥感,2002,23(1):33-35.

[7]杨波,吴德文,赖健清,等.遥感技术在腾冲西南地区地热资源研究预测中的应用[J].国土资源遥感,2003，(2):23-26.

[8]章方,王艳江,杨晓红,等.长春地区地表温度异常遥感分析[J].测绘与空间地理信息,2012,35(9):97-99.

[9]王志才,贾荣光,孙昭民,等.沂沭断裂带安丘-莒县断裂安丘-朱里段几何结构与活动特征[J].地震地质,2005,27(2):212-220.

[10]李小凤.郯庐断裂带[M].北京:地质出版社,2000.

[11]巩贵仁,唐孟武,逯光明,等.遥感技术在临沂地区地热资源勘查中的应用[J].山东国土资源,2008,24(10):17-20.

[12]覃志豪,李文娟,徐斌,等.利用Landsat TM6反演地表温度所需地表辐射率参数的估计方法[J].海洋科学进展,2004,22(B10):129-137.

[13]田国良,柳钦火.热红外遥感(第2版)[M].北京：电子工业出版社，2014.

[14]Chandrasekhar S. Radiative Transfer. New York: Dover Publ Inc,1960.

[15]宫阿都,江樟焰,李京,等.基于Landsat TM图像的北京城市地表温度遥感反演研究[J].遥感信息,2005，(3):18-20.

Application of Landsat ETM+ Remote Sensing Data in Searching Geothermal Field in Linyi City

LIU Xiaocai1，LIU Xiaoyang1, JIANG Yong2, LIU Kun1

(1.No.7 Exploration Institute of Geology and Mineral Resources, Shandong Linyi 276000, China; 2. Linyi Bureau of Land and Resources, Shandong Linyi 276000, China; 3. Land Science and Technology College of China University of Geosciences (Beijing), Beijing 100083, China)

Yishu fault zone is rich in geothermal resources. By using Landsat ETM+ as the main remote sensing data, and inversing the whole surface temperature in Linyi city by using radiative transfer equation. Gray temperature images have been interactively stretched in order to show temperature abnormal areas in 33℃ and the above temperature. Combining with ETM+ true color images and Google images, excluding town, rural residential buildings, river beach, bare rock and other abnormal areas, the already gained abnormal results have been compared with geothermal geologic background images, and analyze it from the angle of regional geology in order to find potential geothermal regions.

Landsat ETM+; temperature inversion; geothermal

2014-07-13；

2014-09-24；编辑：曹丽丽

刘孝才(1983—)，男，山东临沂人，工程师，主要从事地质矿产工作；E-mail:285410775@qq.com

P208

B