玉龙雪山位于我国的云南省丽江市, 呈南北走向, 海拔最高为5500余米, 其温度气候呈现垂直差异, 山顶终年积雪, 区域温差较大。玉龙雪山自然保护区留有我国珍贵的古冰川遗迹, 育有完整自然生态系统, 具有极其重要的科研及旅游价值^[1]。由于其特殊的地理位置, 玉龙雪山对全球气候变化趋势敏感性较强, 因此玉龙雪山的生态环境变化可以很好地反映全球气候变化^[2]。

陆地表面温度是生态环境分析等研究应用的关键指标, 在全球环境变化和生态资源监测等领域具有重要的研究价值^[3]。随着卫星遥感技术的不断发展, 卫星平台的承载力逐渐增大、稳定度不断提高, 为大区域陆地表面温度定量化地反演提供了一种有效的技术手段。目前, 利用卫星热红外遥感影像反演陆地表面温度在地矿资源的探测、城市热岛效应分析、植被病虫害状态监测等领域都取得了重要的研究成果。黄秀华等利用不同时相下热红外遥感数据提取热异常标志, 指出了找油的最佳时机^[4]。张小飞通过提取相关的下垫面类型、地表温度和植被覆盖等信息, 证实城市区域植被覆盖状况可直接影响城市地表温度^[5]。毛克彪提出利用MODIS数据基于辐射传输方程来进行地表温度反演^[6]。周纪构建立了城市热岛容量计算模型, 综合城市热岛强度、足迹面积等多方面信息, 在最大程度上定量地反映和刻画了热岛效应发生的显著程度, 能够客观描述城市热岛对局地气候、人居环境质量的影响^[7]。徐永明利用MOIDS数据提出了改进的温度-植被指数, 在农田区域及农作物生长期内具有很好的适用性和精度, 为有效获取大范围农田气温提供了新的思路^[8]。Brabyn Lars从Landsat 7 ETM +影像提取地表温度, 指出利用Landsat 7 ETM +得出的地表温度可广泛地应用于生态和气象研究^[9]。卢显等基于中国HJ-1B卫星热红外数据验证了2010年4月14日青海玉树M_S7.1地震LST临震地表温度异常的可靠性^[10]。冯海英等人利用MODIS反演的地表温度和土地覆盖产品提出一种评估森林调节温度生态服务价值的新方法, 为正确认识和定量评估森林生态服务价值及进行生态系统管理提供参考^[11]。吴文渊等利用Landsat 8 OLI/TIRS遥感影像反演地表温度, 结果表明地表热环境分布与断裂带等自然因素存在一定相关性^[12]。目前, 基于两线阵、多线阵传感器获取的立体像对重建高精度数字表面模型的研究也取得了很大的进步^[13]。林行刚基于DSM讨论了精确地形对飞行器跟踪制导技术的重要工程应用前景^[14]。王明华采用DSM完成了工程岩体三维地质建模与可视化^[15]。穆超等基于高分辨率遥感影像, 实现了DSM建筑物点的提取, 对建筑物的精确三维重建提供了有效帮助。Ola Friman等采用数字表面模型来估计入射光的不同分量, 进而预测不同光照条件下的被测光谱^[16]。Mathieu Brédifd等提出一种从数字表面模型提取建筑物占地面积的全自动框架, 给出了利用DSM构建3D城市模型的全自动过程^[17]。Lasse Sander等利用SRTM数字地面模型获取到宽滩脊系统的海拔趋势^[18]。凌成星等提出了一种基于GeoEye-1立体像对提取平均树高的方法, 可快速获得研究区大范围森林平均树高^[19]。解金卫等人提出一种高精度DSM估计方法, 能有效提高植被区DSM反演精度^[20]。

热红外遥感能够获取地表温度信息, 可见光遥感具有出色的地面分辨率, 在对地进行精确定位及立体成像方面具有无可比拟的优势。现阶段围绕数字表面模型以及热红外遥感地表温度反演两个领域的应用研究都有了很多重要成果, 但是通过两者的融合应用开展多维遥感信息探测还处于初始阶段, 尚有大量空白空间需要国内外学者进行探索。将星载热红外数据与可见光数据融合使用, 突破单一数据的桎梏, 可进一步拓宽并提升遥感应用的边界与高度, 构建地表三维温度场。相比较单一的热红外遥感应用, 更真实地反应研究区的实际地表情况, 结合温度信息可以实现立体式热红外探测, 通过利用多维信息主被动复合观测, 显著提升目标多维度信息综合监测能力。

本文从星载热红外与可见光立体遥感数据融合应用出发, 以研究玉龙雪山的温度变化为目的, 基于Landsat系列卫星的热红外数据反演二维地表温度, 基于资源三号卫星的三线阵可见光影像构建数字表面模型, 通过两者嵌套完成地表三维温度场构建, 定量分析研究区域的地表三维温度场时空变化。

1 数据与方法 1.1 数据

研究数据选取区域为玉龙雪山区域, 大致位于我国丽江市内100°4'2"E至100°16'30"E、27°3'2"N至27°18'57"N处, 地理位置如图 1所示。研究数据选用陆地卫星(Landsat)系列的热红外影像以及立体测绘卫星资源三号(ZY3)的三线阵影像为基础研究数据。Landsat TM以及Landsat 8 TIRS提供了从1984年至今的热红外数据, 本文选取了从1987年至2018年的时间跨度中, 以两到三年为时间间隔的13个年份中每年的12月热红外影像数据进行处理, 并且通过对比选择该月份中质量较高的日期的影像数据(由于数据质量原因, 在2004年至2008年的数据为四年的时间间隔)。数据在美国地质勘探局网站(USGS)下载得到。资源三号卫星的三线阵影像数据包括前视、后视及全色影像, 轨道行列号Path为33, Row为157, 数据通过自然资源部国土卫星遥感应用中心获取得到。由于较早时期的热红外影像数据无相对应的资源三号数据, 以相近时相下的Landsat ETM+ 全色影像模拟为高分辨率可见光影像对其进行影像配准。以上数据均以Geotiff格式储存, 采用了UTM-WGS84坐标系统进行地图投影。

1.2 方法 1.2.1 普适性单通道算法反演地表温度

普适性单通道算法(generalized single-channel method, GSC)顾及了地表比辐射率和大气辐射参数, 是一种可靠性较高的温度反演算法。其可适用于Landsat 4 TM热红外数据、Landsat 7 ETM+热红外数据和Landsat 8 TIRS数据, 算法公式如下^[21-22]：

(1)

(2)

(3)

式中, T_S为地表温度, 单位为K。T_i为波段通道的亮度温度, 单位为K。L_i为大气顶层卫星传感器所接收到的辐射亮度, 单位为W m^-2sr^-1μm^-1；b_γ=c₂/λ, c₂为普朗克函数中的常量, c₂=1.43877×10⁴μm·K, λ是有效波长。对于Landsat 4 TM第6波段, b_γ为1290K。对于Landsat 5 TM第6波段, b_γ为1256K。对于Landsat 7 ETM+第6波段, b_γ为1277K。对于Landsat8 TIRS第10波段, b_γ为1324K。Ψ₁、Ψ₂及Ψ₃是关于大气水汽含量ω的大气参数, 由下式及表中的计算系数利用大气水汽含量计算得到：

(4)

1.2.2 铅垂线轨迹法制作数字表面模型

铅垂线轨迹法(VLL——Vertical Line Locus)是制作数字表面模型过程中常用的技术策略之一。其基本思想是：假设在物方有一条铅垂线轨迹, 则它在画幅式像片上的投影也是一条直线, 铅垂线与地面交点A在像片上的构像必定位于相应的投影辐射线上, 原理示意图如下^{[13, 23]}。

利用铅垂线轨迹法搜索其相应的像点a₁和a₂, 确定A点高程的基本计算步骤如下：给定地面点的平面坐标(X₀, Y₀)与可能的最低高程Z_min, 高程搜索步距dZ可以按照所要求的高程精度决定。由地面点的平面坐标(X₀, Y₀)与可能的高程：Z_i=Z_min+i·dZ, (i=0, 1, …n), 计算地面点分别在左右像片上的像点坐标(x₁, y₁), (x₂, y₂)。分别以(x₁, y₁)、(x₂, y₂)为中心在各自影像上选取匹配窗口, 计算其匹配测度, 如计算相关系数ρ_i(也可以利用其他测度)。令i=i+1, 重复上述步骤, 得到一组匹配测度{ρ₀, ρ₁, ρ₂, …, ρ_k}, 选取其中的最大值ρ_k, 其对应的高程为Z_k=Z_min+k·dZ, 则认为地面点P的高程值为Z_k。还可以利用ρ_k及其相邻的几个相关系数拟合一条抛物线, 以其极值对应的高程作为地面点的高程, 以进一步提高精度, 或以更小的高程步距dZ, 在一小范围内重复以上过程。

1.2.3 三维温度场构建

地表三维温度场作为星载热红外数据与可见光数据融合产品, 具有几何高精度、信息多维度的优势。在构建地表三维温度场时, 根据应用需求, 在数据选择方面需要考虑热红外数据的时间分辨率和温度分辨率, 例如Landsat系列卫星在理想状态下获取从1984年至今的16天时间分辨率平均误差值在1K左右的地表温度数据, 可广泛应用于各种热红外遥感的研究应用。资源三号卫星三线阵影像具有高精度的几何定位的能力, 采用同轨立体像对生产数字表面模型, 为地表三维温度场提供精确可靠的地面三维坐标数据。

三维温度场的构建基本流程如下：首先对原始可见光影像数据进行预处理, 消除影像的几何与辐射畸变等。采用热红外数据反演地表温度, 得到表征地表温度信息如图 3所示的影像数据。使用可见光数据的前视、后视及正视影像制作如图 3的数字表面模型。最后将热红外温度影像与数字表面模型在统一地理基准的情况下进行套合, 完成热红外与可见光遥感数据地融合, 得到既包含了温度信息又具有高精度三维坐标如图 4所示的地表温度场数据, 技术流程如图 5所示。

3 结果与分析

利用普适性单通道算法进行温度反演, 在获得研究区域的一个时间序列的地表温度反演结果后, 为了更好的研究分析该区域地表三维温度场的时间演变规律及空间分布特征, 需要对地表温度产品进行温度等级的划分。通过对不同温度区间进行分层赋色, 可以直观的表现出该研究区域中不同温度区间的覆盖范围及相对分布情况。本文通过密度分割并基于等间距分割法来划分地表温度等级, 划分低温等级后的温度产品结果如图 6所示：

将以上反演得到的地表温度产品, 与DSM进行融合叠加处理, 获得所有温度产品对应的地表三维温度场。如图 7所示, 不仅可以看出最终生成的地表三维温度场直观明了地展现出玉龙雪山区域在1987年到2018年时间序列中具体地表的温度情况, 同时, 还可以获取到不同地形地物的具体温度信息及三维坐标位置信息。

为了方便对玉龙雪山区域的地表三维温度场进行更直观的定量分析, 根据影像获取时刻的地表温度, 将研究区域简单划分为四个温度区间：低温区(T＜-10℃)、次低温区(-10℃＜T＜0℃)、中温区(0℃＜T＜10℃)以及暖温区(T>10℃)。由于地表三维温度场具有地面及高程的高精度坐标信息, 可以计算得到不同于平面投影面积的真实地表面积, 通过温度属性进行分割可以计算得到不同温区的真实地表面积。对每个温区的地表投影及实际面积进行统计分析对比(表 1、表 2), 并以1987年至2018年中的12月为时间序列绘制不同温区地表面积趋势折线图。

根据上述图表给出的结果, 对玉龙雪山区域从1987—2018年的地表三维温度场的时空变化特征进行研究分析：

(1) 从图 6玉龙雪山区域历年温度平面分布图可以看出, 由于研究区域中部的海拔较高, 温度一致呈现出由中间区域向外递减的趋势, 图 7玉龙雪山区域历年地表三维温度场则更直观的呈现出了历年温度最低的区域一致集中在了海拔最高的山峰顶部, 而在海拔较低及人类活动区域温度逐渐升高。

(2) 由图 8和表 3可以看出, 针对玉龙雪山的时空变化分析, 地表三维温度场相较于二维数据分析优势明显, 对不同温区真实地表面积变化情况更加敏感。其中次低温区位于高程变化显著区域, 真实地表面积与投影面积差异最为明显, 中温区次之。基于地表三维温度场, 充分利用地形高程信息, 可以表现出不同二维数据的真实地表情况。玉龙雪山从河谷到山顶高程差距较大, 对于其垂直气候的分层差异, 表现出不同高程下育有不同的气候带。由于地表三维温度场相比较平面温度场在统计地表面积方面会因高程差异显著更为精确, 因此建立地表三维温度场对于玉龙雪山等山地区域的生态研究分析意义重大。

(3) 在1987—2018年间, 研究区不同温区的地表面积变化呈现出不稳定的变化现象, 整体变化趋于温度升高。由表 2及图 8可知, 所选取的13年中, 1989年玉龙雪山的整体温度分布最低, 低温区(T＜-10℃)所占面积大幅度提高, 之后呈现出逐年递减并在1994年后趋于稳定的状态。次低温区(-10℃＜T＜0℃)从1992年升至最高, 在1989年至1994年相对稳定, 之后表现出逐年下降的趋势。中温区(0℃＜T＜10℃)在1989年降至最低后, 总体表现出先是稳定增加的现象, 在2008年后, 温度分布出现不稳定的状态。相应的, 在1996年, 玉龙雪山的暖温区(T>10℃)地表面积呈现出明显上升, 并在之后的较长的之间序列中变化较小, 但是在2000年后, 不稳定性增强, 年间温度变化较大。从以上现象可以看出, 随着全球温度上升的趋势愈加显著, 玉龙雪山及类似山地或者冰川的自然生态所受影响也愈加明显。同时, 根据相关研究表明, 喜马拉雅山冰川于1990—2015年也呈现出持续退缩状态^[24], 并且其变化会随着高程有所不同。全球变暖造成的影响是整体性的, 同时与高程也密切相关, 本文方法对雪山或冰川相关的生态环境的科学研究有着一定的可行性及优势。

4 结论

本文利用陆地系列卫星和资源三号卫星的影像数据, 通过套合温度反演产品和数字表面模型, 构建了云南省丽江市玉龙雪山区域的地表三维温度场, 研究分析了该区域由1987—2018年的时空变化特征, 并分析不同温区的面积变化。根据本文地表三维温度场与二维数据比较结果, 针对高程起伏明显区域, 地表三维温度场优势明显, 能够更加有效地反映地表温度情况, 对于山地区域的研究具有重要意义。同时, 玉龙雪山区域研究时间序列内的中温区和暖温区的地表面积有逐年增加的趋势, 相较之下中温区的地表面积升幅稍大；低温区和次低温区的地表面积大体呈现减少的趋势；可以发现玉龙雪山地区的整体温度变暖, 验证了全球温度变暖的趋势。