基于太阳辐射效应校正的高原铁路廊道地热异常区识别研究

边 民， 董 庆， 王 栋， 孟德利， 赵文博， 周厚瑀

(1. 中国科学院空天信息创新研究院 数字地球重点实验室， 北京 100094； 2. 中国科学院大学，北京 100049； 3. 中科卫星应用德清研究院 浙江省微波目标特性测量与遥感重点实验室，浙江 湖州 313200； 4. 中铁二院工程集团有限责任公司， 四川 成都 610031)

0 引言

某高原铁路东起成都市，西至拉萨市，是连接我国西南地区的重大工程[1]。其中，新建雅安至林芝段是该高原铁路地质条件最复杂、地形地貌最艰险的地段[2-3]，具有板块活动强烈、生态环境脆弱、地质灾害频发等特点[4]，并且桥梁隧道占比高，极易受到高地温热害的影响[5-6]。高地温热害不仅会降低机械设备效率、影响施工进度，还会恶化作业环境、威胁人员安全，并影响工程后期运营维护[7-8]。因此，准确识别地热异常区对高原铁路建设至关重要。

热红外遥感反演地表温度(land surface temperature，LST)是实现地热异常区识别的一种有效技术手段[9]。但高原铁路沿线地形地貌复杂，地表温度反演会受到太阳辐射效应的影响，高大山坡产生阴影遮挡，使光照面(阳坡)比阴影面(阴坡)受到更多的太阳辐射，地表温度也更高，产生阴阳坡现象[10]，对地热信息的提取带来干扰。因此，去除高原山区地表温度反演中的太阳辐射效应是地热异常识别中亟需开展的工作。

去除太阳辐射的常用方法是通过将存在地形起伏影响的影像像元DN值校正到某一参考平面(通常是水平面)上，抑制太阳辐射效应的影响[11]，考虑到高原复杂的大地形和局部地形多尺度效应影响，太阳辐射效应去除效果欠佳。周桃勇等[12]将地表温度按坡向分成阴坡、阳坡、过渡坡3个子区，然后根据坡向与地表温度的线性拟合关系将其校正到水平坡度上，达到抑制太阳辐射效应的效果，但阳坡、过渡坡和阴坡的划分具有不确定性，且该研究没有对太阳辐射进行定量计算。Hais等[13]将山区植被表层温度表示为海拔和山体阴影的函数，利用线性回归的方法对太阳辐射效应进行削弱，但简单的线性回归模型无法很好地表征出地表温度与各因子之间的关系。已有的太阳辐射校正方法对于太阳辐射没有进行定量的计算，且将太阳辐射对地表温度的影响简化为线性关系，去除太阳辐射效应的效果不理想。随机森林(random forest, RF)作为一种机器学习方法，是目前效果最好的非线性模型之一，可以准确表达地表温度与各因子间复杂的非线性关系[14]；随机森林对多元线性不敏感，无需对变量的正态性和独立性等假设条件进行检验，能够保留样本的最原始信息[15]，并且运算高效、结果准确，可以很好地构建地表温度与各关键因素之间的非线性关系，进而去除地表温度的太阳辐射效应。

除了太阳辐射效应对地表温度反演的影响外，基于热红外遥感识别地热异常区还存在考虑因素不全面、人为主观因素影响大、识别精度较低等缺点。因此，构建一种可靠的地热异常识别模型十分必要。确定性系数(certainty factor，CF)方法是一种二元统计方法，可以用来研究不同因素对某一事件所产生影响的灵敏度[16]；该方法被广泛用于滑坡、泥石流等地质灾害评估[17-18]，具有计算严密、准确性高等优点，可以有效解决输入数据的不确定性和异质性问题[19]。因此，采用确定性系数方法来定量识别和评价地热异常区是可行的。

目前，在使用热红外遥感识别地热异常的研究中，太阳辐射的影响没有得到足够的重视，且仅仅依靠地表温度单一要素，必然会影响对地热异常区的准确识别。本文旨在去除太阳辐射效应引起的地表温度变化，以提高高原铁路廊道地热异常区识别的准确度，然后利用确定性系数方法构建多指标因子模型对地热异常区进行识别和分级，并通过野外温泉点对不同分级的异常区进行定量评价。

1 研究区概况和研究数据

1.1 研究区概况

研究区为高原铁路雅安至林芝段，该区域地处青藏高原东南部(见图1(a))。区域内地质地貌条件复杂，坐落有二郎山、折多山、芒康山、伯舒拉岭、念青唐古拉山等高大山脉，地势陡峭；水系密布，包含大渡河、雅砻江、金沙江、澜沧江、怒江、雅鲁藏布江等水系。受强烈的板块碰撞作用影响，研究区活动断裂(带)发育密集，主要包括鲜水河断裂、理塘断裂、巴塘断裂、金沙江断裂、嘉黎断裂等深大断裂带，具有构造活动强烈、地质灾害多发、温泉地热广布等特点[3]。

(a) 研究区空间范围

(b) Landsat8真彩色影像(2018-11—2019-02)

(c) 研究区DEM

(d) 地磁异常数据铁路线引自文献[3]。图1 研究区与数据概况Fig. 1 Schematic of study area and dataset

1.2 研究数据

本文使用的数据主要包括遥感数据、DEM数据以及其他数据。

遥感数据为白天和夜间的热红外影像。白天热红外数据为2013—2021年冬季覆盖研究区的Landsat8 L1TP级影像共168景(见图1(b))；热红外传感器(thermal infrared sensor，TIRS)数据的空间分辨率已被重采样为30 m，本文进一步对其进行辐射定标、大气校正等预处理。夜间热红外数据来自NASA的空间站生态系统天基热辐射实验(ECOsystem spaceborne thermal radiometer experiment on space station，ECOSTRESS)生成的地表温度和发射率二级数据(ECO2LSTE)，时间范围为2018—2021年冬季；该数据原始空间分辨率为70 m，本文将其重采样至30 m，与Landsat8 TIRS保持一致。

DEM数据来源于ASTER GDEM V3产品(见图1(c))，水平分辨率为30 m。

其他数据包括研究区内温泉点数据(见图1(a))、断裂(带)数据(见图1(a))、地磁异常数据(见图1(d))以及水系数据(见图1(a))。数据的详细情况见表1。

2 研究方法

首先，基于Landsat8热红外数据以及ECO2LSTE数据得到白天和夜间的地表温度；然后，利用随机森林方法去除地表温度反演结果中太阳辐射效应的影响，消除太阳辐射对地热信息提取的干扰；最后，利用多指标因子构建地热异常区识别与评价模型，得到研究区地热异常区分布结果，并结合已知温泉点对结果进行评价。地热异常区识别技术流程图见图2。

图2 地热异常区识别技术流程图Fig. 2 Flowchart of geothermal anomaly identification

2.1 地表温度反演

地表温度的反演算法主要有辐射传输方程法、通用单通道算法、单窗算法、分裂窗算法、温度发射率分离法等[20]。由于Landsat8卫星第11波段的定标参数不稳定，本文使用单窗算法[21]反演地表温度。

单窗算法基于辐射传输方程推导得到，只需要已知地表发射率、大气透过率和平均温度3个参数就可以反演得到地表温度：

ts=[a(1-C-D)+b(1-C-D)+C+D)tb-Dta]/C。

(1)

式中：ts为卫星反演得到的地表温度，℃;ta为大气平均温度，℃；tb为亮度温度，℃；a和b为回归系数,分别取值-62.72和0.433 9;C和D为中间变量，可通过大气透过率和地表发射率计算得到[21]。

2.2 太阳辐射效应去除

本文采用随机森林方法和白天/夜间的多时相地温求均值的方法，对地表温度的太阳辐射效应进行抑制，突出地热造成的地表温度异常。

2.2.1 随机森林方法

地表温度与地形、植被条件、水分条件和太阳辐射密切相关。本文采用随机森林方法构建模型，以准确表达地表温度与关键要素之间的非线性关系，模型的自变量为海拔、坡度、坡向、累积太阳辐射、NDVI(归一化植被指数)、NDSI(归一化积雪指数)、NDWI(归一化水体指数)以及反照率(见图3)，因变量为地表温度。累积太阳辐射采用Han等[22]的方法计算得到，反照率采用Traversa等[23]的方法计算得到。利用模型模拟地表温度与原始反演地表温度进行对比可以验证模型模拟的准确度，基于随机森林模型对各因子与地表温度非线性关系的准确模拟，进而将模型中与太阳辐射效应相关的输入因子(如坡度、坡向和累积太阳辐射等)设为参考值，而其余输入因子保持不变，得到去除太阳辐射效应的地表温度。

(a) 海拔

(b) 坡度

(c) 坡向

(d) 累积太阳辐射

(e) NDVI

(f) NDSI

(g) NDWI

(h) 反照率图3 随机森林模型自变量Fig. 3 Independent variables of random forest model

1)通过随机森林回归算法建立地表温度与各影响要素之间的非线性模型：

tRF=F(ELV,SLP,ASP,CSR,NDWI,NDVI,NDSI,ALB)。

(2)

式中：tRF为使用随机森林回归算法模拟的地表温度，℃；F为非线性函数模型； ELV为海拔高程，m； SLP为地表坡度； ASP为地表坡向； CSR为累积太阳辐射，MJ/m2； ALB为地表反照率。

2)随机森林方法模拟的地表温度与原始地表温度之间会存在一定的误差，可以表示为：

Δ=tO-tRF。

(3)

式中Δ为原始地表温度值tO与随机森林模拟温度值tRF的差值，可以认为是随机森林模型的模拟误差。

3)结合式(2)和式(3)，可以得到：

tO=F(ELV,SLP,ASP,CSR,NDWI,NDVI,NDSI,

ALB)+Δ。

(4)

地表温度值可以看作是多种关键要素的函数，其中，包含由太阳辐射差异引起的变化值和非太阳辐射相关值。本文通过随机森林模型去除地表温度的太阳辐射效应，得到校正后的地表温度。

2.2.2 日夜多时相地温求均值

利用日夜多时相地表温度联合探测地热异常，可以更好地去除太阳辐射效应对地热异常识别的干扰。夜间地表温度虽不受太阳辐射的直接影响，但河谷效应显著；河谷效应也是由太阳辐射效应引起的，因为水的比热容大，且白天水体受到太阳辐射蒸发，河谷内空气湿热，夜间由于河谷地带地形相对封闭，加剧了地表升温，也会对地热信息的提取带来干扰。由于不同地物比热容及赋存结构的差异，自然地物在白天/夜间的相对热/冷异常对地热信息提取的干扰也必须考虑在内。因此，本文采用日夜多源多时相热红外遥感数据相加求平均值，通过白天和夜间温度联合分析，消除太阳辐射效应和自然地物间的相对热/冷异常，更好地突出地热所导致的地温高异常，提升对地热异常区的识别能力。

2.3 地热异常区识别

地热高温异常的产生是多种因素共同作用造成的，需要综合考虑不同因子对地热分布的影响，并确定指标量化方法，选取适当的数学模型计算得到地热异常的易发程度。本文采用确定性系数(certainty factor，CF)模型来识别和评价地热异常区。

2.3.1 确定性系数值计算

确定性系数是一个概率函数[24]，主要思想是地热的易发性可依据已经存在的地热信息与确定为影响因子的数据集之间的统计关系来确定[25]。地表出露的温泉是地热异常的有效表征，因此，可将温泉视为已知的地热信息，参与地热异常识别模型的构建。由于影响地热存在的各指标因子量纲不同，不能直接进行数学计算，而确定性系数可将各因子同区间定量化，有效解决不同数据量纲的差异性问题[21]。确定性系数计算公式为：

(5)

式中：KCF为确定性系数；Pa为影响因子分级图层内存在地热的条件概率，此处为分级内温泉点数目与分级内栅格数目之比；Ps为地热发生的先验概率，此处为研究区有温泉点的评价单元数与研究区评价单元总数之比。

利用确定性系数模型识别和评价地热异常区，需要对选取的指标因子进行分级，并计算各分级区间内的确定性系数大小。确定性系数的取值范围是[-1,1]，值越大代表单元内存在地热的确定性越高。当结果为正时，表示有利于地热异常的存在；当结果为负时，表示不利于地热异常的存在；当结果接近于0时，无法确定评价单元是否存在地热异常。

2.3.2 各指标因子权重计算

从式(5)可以看出，各指标因子的确定性系数值与该因子对地热存在与否的贡献值密切相关，因此，可以根据CF值计算每个因子的相对贡献大小。

(6)

式中：w为某指标因子的相对贡献值；Zi为某因子的CF分段贡献值；i为CF的分段级别，i=1，2，…，n。

将式(6)计算的各指标因子相对贡献值归一化后，即可得到该因子的权重。

3 研究结果与分析

3.1 地表温度反演结果

研究区多年期冬季平均地表温度分布如图4所示。由于本次研究针对的是高地温异常，故将白天和夜间地表温度的下限值分别设定为-20 ℃和-30 ℃；同时，将白天和夜间地表温度的上限值分别设定为40 ℃和20 ℃，温度区间统一，便于对比太阳辐射效应校正前后地表温度的差异。可以看出，Landsat8反演得到的白天地表温度受到地形背景和太阳辐射效应的显著影响，地表温度高值区大多分布于地形起伏较大区域的山坡阳面、河谷地带以及裸露地表区域；低值区主要是积雪、山体阴影以及植被密集区域。ECO2LSTE夜间地表温度的河谷效应显著；高值区主要是水体和河谷地区，大致包括大渡河、雅砻江、金沙江、澜沧江、怒江、雅鲁藏布江等大江大河的河谷地带；低值区主要是裸地、积雪覆盖地区以及高海拔的山区。

黑框区域为阴阳坡效应典型区。(a) Landsat8白天地表温度

黑框区域为河谷效应典型区。(b) ECO2LSTE夜间地表温度图4 研究区多年期平均地表温度Fig. 4 Multi-year average LST of study area

3.2 地表温度太阳辐射效应去除

在随机森林模型中，ntree是对模型效果影响程度最大的参数，需要择优选择。ntree值和模型预测得分如图5所示。从图中可以看出，当ntree值为0～25时，随着ntree取值增大，模型的预测得分也会增加，且增长速率较快；当ntree值为25～100时，模型的预测得分增长速率逐渐放缓；当ntree值达到100后，模型预测得分趋于稳定。因此，本文ntree取值为100。

图5 ntree值和模型预测得分Fig. 5 Plot of ntree number and model prediction score

随机森林方法去除太阳辐射效应的前提是基于对地表温度的精确模拟。图6(a)和图6(b)分别示出白天和夜间模型模拟LST的空间分布，并将局部的模型模拟LST与原始反演LST进行对比验证。从图中可以看出，模型模拟LST与原始反演LST的一致性程度较好，LST高值区和低值区高度吻合。对局部放大图进行对比分析，模型模拟的LST可以较好地反映LST的层次过渡且能够很好地表征阴阳坡效应、河谷效应等现象，表明随机森林模型模拟的精度较好。

(a) Landsat8白天模拟地表温度

为进一步定量分析原始反演LST与模型模拟LST之间的一致性，本文采用决定系数(R2)、均方根误差(RMSD)、平均绝对误差(MAD)作为评价指标，来评价随机森林模型的模拟效果。原始反演LST与模型模拟LST比较如图7所示。训练数据的比较结果显示，影像的R2为0.97，RMSD为0.93 ℃，MAD为0.66 ℃；测试数据的比较结果显示，影像的R2为0.92，RMSD为1.95 ℃，MAD为1.76 ℃。上述统计参数表明，模型具有较高的模拟精度，随机森林方法可以精确解译地表温度空间分布特征。

(a) 训练数据

(b) 测试数据散点的颜色代表了该位置点的密度。图7 原始反演LST与模型模拟LST比较Fig. 7 Comparison result of original inversion LST and model simulation LST

太阳辐射效应去除后的地表温度分布如图8所示。研究区由于地形背景和太阳辐射差异造成阴阳坡效应、河谷效应明显，地表温度分布变化剧烈，而去除太阳辐射效应的地表温度与原始地表温度相比变化范围明显减小，LST空间异质性明显降低，太阳辐射得到了较明显的消减。分别选取研究区内阴阳坡效应和河谷效应典型区域进行分析。1)样区1位于雅江县城东部，区域内地形起伏明显，且山地呈东西走向，具有显著的阴阳坡地理分异。对比校正前后的LST分布图可以看出，校正前的LST空间差异十分显著，太阳辐射效应尤其突出；而校正后的LST分布较为均匀，过渡更为平滑，空间异质性明显降低，阴坡阳坡的地表温度差大幅减小。2)样区2位于雅江县城所在河谷地区，该区域河谷深切，地表温度具有典型的河谷效应。对于夜间LST，校正前河谷地区地表温度显著高于周边地区，而校正后河谷地区地表温度的高值区得到了有效抑制。由此可见，基于随机森林的LST地形效应去除方法效果显著。

(a) Landsat8白天校正后地表温度

(b) ECO2LSTE夜间校正后地表温度图8 太阳辐射效应校正LSTFig. 8 LST after solar radiation effect correction

经过太阳辐射效应去除后的日夜多时相地表温度平均值可以作为后续确定性系数模型的输入因子，用于地热异常区的识别。

3.3 地热异常区识别与评价

3.3.1 确定性系数计算结果

地热的产生会受到多种因素的影响，基于地热异常的形成机制和分布规律，结合研究区的地形、地质和实地调查资料，本文选取去除太阳辐射效应的地表温度、断裂密度、到水系距离和地磁异常4个评价指标作为模型输入因子。地热蕴含着丰富的能量，其引发的地球表面温度异常，可利用热红外遥感反演来获取，因此地表温度是地热存在的重要指示因子[26]；地下水可以从断裂形成的渗流通道到达地壳深处，加热后产生地热活动[27]，基于断裂(带)数据可以得到断裂密度作为模型输入因子；在一定深度范围内，地表水和浅层地下水沿着构造裂缝渗流到深部，在下渗过程中被围岩加热而形成地热[28]，因此水系可为地热的形成提供水源，采用到水系距离作为地热分布的影响因素是可行的；地磁异常可以用来表征构造应力变化明显区域以及地下热水活动区域，因此可将地磁异常作为输入因子[29]。将选取的因子重分类为7个等级(见图9)，并计算各因子的确定性系数值(见表2)。

(a) 地表温度

(b) 断裂密度

(c) 到水系距离

(d) 地磁异常图9 确定性系数模型输入因子Fig. 9 Certainty factor model input factor

表2 确定性系数计算结果Table 2 Calculation results of certainty factor

由式(6)计算得到各因子的相对贡献值并归一化后得到地表温度、断裂密度、到水系距离和地磁异常的权重分别为0.35、0.28、0.24、0.13。将各因子的确定性系数值加权叠加获得CF模型最终计算结果，如果模型值低于0，则说明与地热高温现象无关，因此，将第1个阈值设置为0，模型大于0的值按照自然断点法分为3类，值越高代表存在地热异常可能性越高。最终将结果划分为高异常区、中异常区、低异常区和无异常区4个区域，如图10(a)所示。高异常区和中异常区主要包括康定县北部和西部、道孚县中部和北部、理塘县中部、贡觉县和察雅县西部、八宿县西部以及林芝县东北部地区，这些地区受到鲜水河断裂带、甘孜—理塘断裂带、金沙江断裂带、怒江断裂带、澜沧江断裂带以及嘉黎断裂带的控制，且多位于河谷地区，具有地热产生的条件。

文献[30]对高原铁路廊道地热异常区进行了划定，对地热异常区的空间分布特征和铁路沿线隧道热害进行了定性识别。将其结果(见图10(b))与本文得到的地热异常区进行比较和评价，可以看出本文与前者对地热异常区的识别基本吻合，但本文的异常区分级更符合实际情况。从异常区分布来看，前者对异常区的识别范围过大，包含了很多非异常区，如丹巴县西部、八宿县南部、波密县东部以及昌都和察雅交界地区。造成这一现象的主要原因是文献[30]按照不同控热构造及温泉分布综合解译划分异常区，以定性识别为主，识别结果空间分辨率较低，结果的可靠性依赖于野外实测资料，具有一定的主观性和不确定性。而本文中地表温度去除了太阳辐射的影响，采用多指标因子进行定量化识别，划定的异常区范围更加精细，识别结果可靠性更高。

(a) 本文所得地热异常区分级

(b) 文献[30]所得地热异常区分级图10 地热异常分级图Fig. 10 Geothermal anomaly classification map

3.3.2 铁路沿线热害分布

相比于整个研究区，铁路沿线两侧是需要重点识别高温热害的区域，因此，本研究分别根据高原铁路沿线左右20 km和30 km作为缓冲区(内层20 km为核心区，外层30 km为次级区)与地热异常区进行叠加分析(见图11)，并对不同等级地热异常区交通廊道长度进行分析和统计(见表3)。可以看出，高原铁路主要途径泸定—康定、理塘盆地及周边、巴塘、贡觉、昌都—察雅、八宿、波密—鲁朗7个地热异常区。

红色虚线圈定区域为地热异常区。图11 铁路沿线缓冲区地热分布Fig. 11 Geothermal distribution in buffer zone along railway

表3 高原铁路廊道地热异常区统计Table 3 Statistics of geothermal anomaly areas along plateau railway

3.3.3 地热异常提取结果评价

频率比可用于评估模型识别结果的有效性和可靠性[17]。该方法是将识别结果每一分级内的温泉占比除以该级别内栅格占比得到频率比，可说明模型的识别结果是否准确。从频率比计算结果(见表4)可以看出： 1)高异常区栅格数仅占研究区总栅格数的4.40%，却包含有66个温泉点，占温泉点总数的33.50%，频率比达到7.61；中异常区栅格数占研究区总栅格数的8.35%，涵盖51个温泉点，占温泉点总数的25.89%，频率比为3.10；高异常区和中异常区的频率比均高于1。2)无异常区的栅格占比达74.00%，而温泉点占比仅为16.75%，频率比为0.23，低于1。3)从低异常区到高异常区，随着异常程度的升高，频率比逐渐增大(见图12)，这说明模型结果是合理可靠的，能够准确识别地热异常区。

表4 频率比计算结果Table 4 Frequency ratio results

图12 不同地热异常分级频率比变化Fig. 12 Frequency ratio variation of different geothermal anomaly areas

4 结论与讨论

山区太阳辐射效应的影响给热红外识别地热异常造成困难，针对这一问题，本文以高原铁路廊道为研究区，利用白天和夜间地表温度数据，基于随机森林方法去除地表温度的太阳辐射效应，结合断裂数据、温泉点与水系数据，利用确定性系数模型识别地热异常区并对其结果进行了定量评价，得到以下结论：

1)随机森林方法可以较好地模拟Landsat8反演的地表温度，并且能够有效消除太阳辐射效应。校正之后的地表温度空间异质性显著降低，太阳辐射影响造成的LST异常得到了有效抑制，特别是在高程起伏较大、阴阳坡和坡度变化大的山区，太阳辐射效应去除效果更佳。

2)太阳辐射效应去除后的地表温度可以更为有效地识别和提取地热异常区。选取去除太阳辐射效应的地表温度、断裂密度、到水系距离和地磁异常作为输入因子，利用确定性系数模型计算得到地热异常区，并将结果分为高异常区、中异常区、低异常区和无异常区4个等级。地热异常区评价结果与已知温泉点分布较为一致，表明模型结果合理可靠。

3)高原铁路主要途径7个地热异常区，地热异常区的统计结果表明高原铁路经过高异常区的廊道长度为72 km，占比7.1%；中异常区的廊道长度为125 km，占比12.3%。中异常区和高异常区主要分布于泸定—康定、理塘盆地及周边、巴塘、贡觉、昌都—察雅、八宿、波密—鲁朗，这些区域水系密布、多深切河谷，断裂活动性强，断裂与断层破碎带作为导水构造，为深部热源的循环运移提供了有利条件，地表温度较高，是高原铁路施工过程中需要重点防范高温热害的区段。

本文研究了如何抑制地表温度反演中太阳辐射效应影响，并对高原铁路沿线的地热异常区进行了定量识别和分级评价，初步实现了高原铁路廊道的地热异常区定量化和精细化识别，可为山区高地温热害等工程地质风险的早期识别提供科学依据。但是，仍然存在不足之处： 一方面，除了太阳辐射效应外，地表覆被和地貌对地热异常的识别也会有重要的影响；另一方面，考虑到夜间温度数据较少、分辨率不高，且山区河谷地表温度成因复杂，河谷效应去除仍需改进和提升。因此，在接下来的研究中将构建太阳辐射影响的连续地表温度模型，重建同期的夜间温度数据，以更好地去除河谷效应和地表覆被的影响，进一步提升地热异常区的早期识别可靠性和准确性。

致谢

感谢中国-巴基斯坦地球科学研究中心提供研究支持。