基于多源RS数据的多尺度青藏高原冻土监测

黄博文，高瑞翔，陈一仰，范增辉，刘修国

(中国地质大学(武汉)地理与信息工程学院，湖北 武汉430074)

冻土是指土壤温度在零摄氏度及以下，并含有冰的各种岩石和土壤，一般可分为短时冻土(数小时/数日以至半月)、季节冻土(半月至数月)和多年冻土(又称永久冻土，指的是持续2年或2年以上的 冻结不融的土层)。在全球变暖的影响下，冻土已经成为全球气候变化的重要观测对象之一，通过对冻土的观测可以一定程度上反映全球变暖的严重程度[1]。

青藏高原作为世界上中低纬度海拔最高、面积最大的多年冻土区域，其冻土变化会对中国东部乃 至东亚气候的形成、变化和发展造成重大的影响。由于青藏高原海拔高、地域广、环境恶劣，通过地面调查研究冻土环境的工作方法难以大范围地应用。许多基于温度数据的冻土分布模型往往需要依靠地表温度观测站点获取地表温度数据，具有空间局限性明显、监测密度低、人力物力成本较高等特点，也难以在区域尺度上应用。遥感(RS)技术能够在一定的条件下进行大范围的数据采集，较好地解决了由于环境险恶造成的交通不便、测量困难等难题，弥补了传统研究技术上的不足。本文针对青藏高原冻土的特点，提出一种基于多源RS数据融合技术反演地表温度，并结合相关的算法，实现冻土空间分布提取和分析的方法。该方法可以消除极端天气对地表温度反演结果的不利影响，实现冻土环境的精细化监测。

1 研究区概况与RS数据源

1.1 研究区概况

本文研究区为青藏高原地区。青藏高原位于北纬25°～40°E、东经74°～104°E之间，面积约250万km2，平均海拔为4 500 m。主要地形地貌为山地，其中包括东西与南北走向的两种山脉，其间还分布大量的湖泊与河流。该地区气温常年较低，部分地区年积温低于0℃，为冻土形成创造了先决的条件。 青藏高原多年冻土区面积约为157万km2,主要以亚稳定型(30.4%)、过渡型(22.1%)和不稳定型(22.6%)多年冻土为主，极稳定型和稳定型多年冻土的占比较少。此外，选择黑山实验区为重点研究区。它位于冻土北界，库塞湖东北方向，在柴达木盆地东北边缘，深藏于祁连山腹地。该实验区毗邻黑海，地形多样,地表覆盖类型众多,包括植被、冰川、山地、湖泊等，有利于检验RS监测中多源数据快速融合并提取冻土的方法。

1.2 RS数据源的选择及其特征

目前国内外提出了许多基于温度数据的冻土分布模型如冻结指数模型[2]，还有利用SAR数据进行冻土地表覆盖提取[3]等。随着遥感科学的发展，通过遥感影像反演地表温度的技术已经较为成熟，中分辨率成像光谱仪(MODIS)温度产品数据已在冻土监测研究中得到了应用[4]。赵全宁等[5]通过线性回归分析指出，在未来气候变暖背景下，在玉树地区应用温度因子估算平均最大冻土深度的变化具有较高的可信度；秦艳慧等[6]对ERA-Interim进行了重新拟合改良,并使用其通过拟合得到了我国青藏高原地区的冻土分布图，与实际情况对比其拟合效果较好。此外，近年来的相关研究充分说明了使用遥感相关的温度产品反演青藏高原地区的冻土分布是可行的。

本文选择了MODIS的8天合成地表温度数据(MOD11A2)作为大尺度温度数据，该数据时间分辨率高，涵盖白天和夜晚地表温度信息，采用的是8天内晴好天气状况下的地表温度数据，并根据地表覆盖类型数据进行平均计算而得到，能够较好地反映地表在一段时间内的温度分布情况,其空间分辨率为1 km。依据该数据可以计算出相应的年平均地温和年积温，用于冻土分布信息的自动化提取。小尺度影像选择Landsat8 影像的2级数据，该等级数据经过辐射校正和几何校正，空间分辨率为30 m，重返周期约为16 d,图像的精细度较高，可以利用红外波段数据反演得到地表温度结果。此外，水汽数据来源于MOD05-L2系列产品,其反映的是晴空区域陆地上空的云层上方的大气柱水蒸气量，通过该数据以及Landsat8的影像数据即可根据算法得到地表温度。

地表覆盖类型数据采用:①清华大学制作的全球地表覆盖类型数据，该数据由清华大学宫鹏教授团队制作，包括2017年全球0.000 25度(对应约30 m)土地覆盖空间分布，包含农田、森林、湿地、水体 等10类地表覆盖类型；②美国NASA制作的全球1 km地表覆盖类型分类。这两种产品数据精度较高、分辨率尚可，经过裁剪处理后可以作为中间变量直接使用。

2 地表温度数据时空融合方法

2.1 数据预处理

2.1.1 MODIS数据预处理

青藏高原区域涉及6景MODIS数据，经拼接裁剪后方可得到研究区的完整影像。全年的MOD11A2数据共计276景影像。本研究的时间跨度为2000年到2018年,计19年，共需要处理5 244景影像。使用IDL语言调用MCTK插件对MODIS数据进行重投影、裁剪、镶嵌等批量化预处理，较快地得到了MOD11温度产品数据集。

由于大气影响和云雾遮挡，MOD11A2数据中存在着异常值和无效值等随机误差，可依据拉依达准则[7]进行误差消除。另外，由于地表温度受太阳直射的影响大,仅使用日温度数据将会导致反演的年平均地温过高，因此需综合使用MOD11A2的日夜温度数据，采用日夜温度平均的方式反演可信的地表温度数据[4]。

2.1.2 Landsat影像数据预处理

Landsat8数据的预处理包括：裁剪、拼接、辐射定标和去云处理。使用ENVI软件中的内嵌工具即可完成Landsat8数据的预处理。对于有云的部分，利用Landsat8数据提供的云层数据，对影像进行掩膜处理，可以获取去云后的 Landsat8数据。同时，可根据该数据提供的反射率去除高海拔的冰川，可避免冰川被误分类而影响结果精度。

2.1.3 地表覆盖类型数据预处理

对于下载的地表覆盖类型数据文件，进行拼接、裁剪处理，可获得研究区内的地物分类。该数据的主要作用是确定不同地表覆盖类型下的土壤导热系数，以提高地表温度反演与冻土类型确定的计算精度。

2.2 Landsat与MODIS地表温度数据的融合方法

2.2.1 Landsat地表温度反演

传统的单通道算法[8]是通过估计大气参数修正单一通道中的大气效应来反演地表温度。Practical Single-Channel(PSC)算法[9]在单通道算法的基础上，针对不同传感器计算固定参数，降低了单窗算法中由于普朗克函数线性化和大气校正带来的误差。该算法分两步进行地表温度TS估计：首先通过消除地表发射率和大气效应引起的误差，计算地表黑体辐射亮温B(TS)；然后根据普朗克函数的反函数，计算地表温度TS。具体实现过程为通过公式(1)、(2)计算地表比辐射率，再通过公式(3)计算地表黑体辐射亮温，将计算得到的地表黑体辐射亮温代入公式(4)即可计算地表温度结果。详见公式(1)～(4)：

(1)

(2)

(3)

(4)

式中：PV表示植被覆盖度，其中NDVI表示归一化差异植被指数；NDVIS表示裸地的归一化差异植被指数(取值为0.2)；NDVIV表示纯植被的归一化差异植被指数(取值为0.55);w表示水汽数据，来源于MOD05数据;ε表示地表比辐射率；εS表示裸地的比辐射率(取值为0.966 8)；εV表示纯植被的比辐射率(取值为0.986 3)；Lsen表示星上辐射亮度，为采用Landsat8的热红外波段辐射定标之后的结果;参数a0～a7均为与传感器有关的固定参数，本文使用Landsat配套参数，具体见表1；λ表示中心波长(nm)；C1= 1.191 04×108W·μm4/(m2·sr),C2=1.438 77×104μm·K。

表1 Landsat配套参数表Table 1 Parameter-value of Landsat

Landsat遥感影像在研究区的成像时间为中午时分，依照此方法反演出的地表温度要比MOD11A2的日夜平均温度高，使得最终的年平均地温的计算结果较使用MOD11A2数据高出许多，是一个极大的不可忽略的误差因素。不同RS数据源获取的对应时间段内的地面温度相互关联，可以通过拟合分析技术研究其相关性[10]。本文采用了多项式拟合技术，将Landsat遥感影像数据反演出的地表温度与对应的MOD11A2地表温度数据进行拟合,MOD11A2数据的时相与Landsat遥感影像相一致，使得拟合结果更为可靠。采用三次多项式拟合，拟合结果即为地面当日平均温度，将拟合后的夏季与冬季温度作为年最高与最低温度计算年平均地温。

2.2.2 冻融指数计算

冻融指数是冻土判识与类型区划的依据，本文采用的计算方法与文献[11]相同，其关键参数是确定累计正温DDT和累计负温DDF。基于地温全年变化符合余弦曲线规律的假定，DDT和DDF的计算公式如下：

(6)

(7)

(8)

(9)

(10)

式中：Th代表全年最高月温(℃)；Tl代表全年最低月温(℃)。

根据上述公式，只需使用全年极值温度进行计算即可，具体计算示意图见图1。此种处理方式的优点在于数据处理量小，一年仅需要两个极值月份的影像即可以计算冻融指数。根据冻土北界黑山实验区内五道梁站点的地表温度观测数据[12]，区内最低温度所在月份为1月或12月，最高温度所在月份为7月或8月。在需要精细化冻土分布的区域，本文利用该计算方法对拟合处理后的Landsat温度数据进行冻融指数的估算。

图1 地表温度拟合曲线示意图Fig.1 Sketch of surface temperature fitting curve 图中：A为拟合曲线的振幅;Tp为拟合曲线的平均值所在处，可以参考图中所标注信息；DDT为累计正温(℃)；DDF为累计负温(℃)。

由于全球气候变化的影响，气候波动和极端气候出现的次数有所增加，使得年内温度变化不再符合简单的简谐波动假设，会造成积温计算的不准确。 为此，本文提出一种针对MODIS全年温度产品数据集的冻融指数计算方法，以0℃为界限，将一年的温度划分为正温和负温，并分别累加计算正积温和负积温,从而得到准确的冻结指数和融化指数[13]。该方法从根本上消除了气候波动对年内积温计算准确性的影响，实验表明该方法计算得到的年平均地温的精度更高，有利于后续的冻土分类与制图。

2.3 基于温度的冻土类型识别与分类

以多源RS数据为基础，本文利用TTOP(Temperature at the Top Of Permafrost)模型反演活动层顶板的温度(TTOP)，其计算公式如下[11]：

(11)

式中：rk为冻融指数；p为土壤导热系数[W/(m·K)]。

根据公式(11)，通过输入土壤导热系数和冻融指数，即可反演得到活动层顶板的温度。利用反演出的活动层顶板温度，根据程国栋院士提出的热力学分布方法[14]，可以利用活动层顶板温度实现冻土分类与制图。

本文选取MODIS温度产品2001年到2018年共18 年的数据，根据TTOP模型反演得到每年活动层顶板的平均温度，对所有数据在时间域上进行滤波处理,以降低异常年份对冻土反演结果的影响。将数据分为两期，在每期结果内对反演得到的活动层顶板温度进行平均，据9年平均结果得到最终的TTOP作为分类依据，见表2。

表2 冻土类型划分标准Table 2 Permafrost classification standards

3 地表温度数据融合和冻土类型RS识别

3.1 Landsat与MODIS地温数据融合结果

利用上述Landsat与MODIS数据的融合方法，可得到2000—2018年青藏高原Landsat8地温反演结果。本文任选2013年和2018年青藏高原Landsat8地表温度的反演结果进行分析，见图2和图3。

图2 2013年青藏高原Landsat8地表温度的反演结果Fig.2 Inversion result of surface temperature of Qinghai- Tibet Plateau in 2013 based on Landsat8

图3 2018年青藏高原Landsat8地表温度的反演结果Fig.3 Inversion result of surface temperature of Qinghai- Tibet Plateau in 2018 based on Landsat8

3.2 基于地表温度的冻土类型RS识别结果

多年年平均活动层顶板温度可由活动层顶板温度TTOP数据经平均化处理后得到，进而反演得到青藏高原的冻土分布图，见图4 和图5。

图4 2001—2009年青藏高原冻土分布图Fig.4 Distribution map of permafrost of Qinghai- Tibet Plateau in 2001—2009

图5 2010—2018年青藏高原冻土分布图Fig.5 Distribution map of permafrost of Qinghai- Tibet Plateau in 2010—2018

3.3 精度评价

3.3.1 与历史数据的一致性分析

将2001—2009年青藏高原冻土分布图(见图4)和2010—2018年青藏高原冻土分布图(见图5)分别与1∶300万青藏高原冻土分布图(见图6)矢量数据和青藏高原历史冻土分布图(见图7)进行宏观上的叠加比对，其结果见图8和图9。结果发现：1∶300万青藏高原冻土的分布范围与图4和图5基本一致，说明在宏观上采用本文的冻土环境RS信息提取的技术方法是可信的[15]；在多年冻土的集中分布区域(如图6中箭头所示)，青藏高原历史冻土分布(见图7)与图4、图5基本吻合;2001—2009年青藏高原冻土分布与历史冻土分布的相似性达90.42%(见图8)，2010—2018年青藏高原冻土分布与历史冻土分布的相似性达90.13%(见图9)。由图10和图11可见，利用Landsat所反演出的2013年青藏高原冻土分布的相似性约为76.24%(见图10)，2018年冻土分布的相似性约为74.57%(见图11)。

图6 1∶300万青藏高原冻土分布图Fig.6 Distribution map of permafrost of Qinghai- Tibet Plateau in the scale of 1∶3 000 000

图7 青藏高原历史冻土分布图Fig.7 Comparison between the distribution of permafrost of Qinghai-Tibet Plateau in 2010—2018 and historical permafrost

图8 2001—2009年青藏高原冻土分布图与历史 冻土分布图的叠加对比Fig.8 Distribution map of historical permafrost of Qinghai-Tibet Plateau

图9 青藏高原2010—2018年冻土分布图与历史 冻土分布图的叠加对比Fig.9 Comparison between the distribution of permafrost of Qinghai-Tibet Plateau in 2001—2009 and historical permafrost

图10 利用Landsat8反演的2013年青藏高原冻土分布图Fig.10 Distribution map of permafrost of Qinghai- Tibet Plateau in 2013 based on Landsat8

图11 利用Landsat8反演的2018年青藏高原冻土分布图Fig.11 Distribution map of permafrost of Qinghai- Tibet Plateau in 2018 based on Landsat8

3.3.2 Landsat的年平均地温拟合的准确度

以青藏高原2013年和2018年度为例，将Landsat反演得到的年平均地温数据与MODIS标准温度进行了对比，其差值结果见图12和图13。

图12 利用Landsat反演的2013年青藏高原年平均地 温与MODIS标准温度的差值分布图Fig.12 Distribution map of the difference between the mean annual ground temperature and MODIS standard temperature of Qinghai- Tibet Plateau in 2013 based on Landsat8

图13 利用Landsat反演的2018年青藏高原年平均地 温与MODIS标准温度的差值分布图Fig.13 Distribution map of the difference between the mean annual ground temperature and MODIS standard temperature of Qinghai- Tibet Plateau in 2018 based on Landsat8

由图12和图13可见，利用Landsat反演得到的2013年青藏高原年平均地温与MODIS标准温度的平均误差为0.419，方差为3.01；2018年相应的数据与MODIS标准温度的平均误差为0.619，方差为2.34，说明本文采用的融合技术方法获取的年平均地温数据的准确度较高。

4 结论与建议

本文针对青藏高原冻土的特点，提出了一种利用多源遥感数据进行地表温度提取的反演方法。该技术方法通过使用MODIS 8天温度合成数据与Landsat8 反演出的地表温度进行拟合，最终得出地表的年平均温度；再根据TTOP模型计算冻融指数，获得以Landsat8为基准的冻土分布，其空间精细程度与准确度相较于以往的大尺度制图有了较大程度的提高，研究成果可为相关领域的应用提供基础数据，技术方法也可为类似研究提供借鉴。

由于地表温度反演受到多方面因素的影响，其准确度仍然存在一定的偏差，因此还需要进一步考虑地表温度反演分析的影响因素及其相关误差的来源，在此基础上提高遥感数据的提取精度。