2017年中巴经济走廊冻土分布数据集

艾鸣浩，张耀南＊，康建芳，冯克庭，田德宇

1. 中国科学院西北生态环境研究院，科学大数据中心，兰州 730070

2. 国家特殊环境、特殊功能观测研究台站共享服务平台，平台服务中心，兰州730070

数据库（集）基本信息简介

引 言

冻土是北半球冰冻圈分布最广泛的因子，占陆地面积的56%，其中多年冻土约占24%[1]。多年冻土具有热稳定性差、水热活动强烈、厚层地下冰和高含冰量冻土所占比重大等特点，在全球气候变暖背景下，多年冻土对工程活动的热扰动极为敏感。在多年冻土区工程建设中，施工期不可避免地会改变包括地表覆盖、地形因子及浅层地表能量平衡等局地因素，加上开通运营后动载荷对冻融循环水热环境的改变，直接导致冻土温度升高和上限抬升，进而增加了冻胀、融沉和冻拔等冻土灾害风险，对多年冻土区工程施工和安全运营构成严重威胁，也对区域水文地质条件和生态环境演化产生巨大影响[2-4]。

中巴经济走廊作为“一带一路”倡议的重要组成部分，是中国与中亚、西亚、南亚地区的交通要道，对中国和巴基斯坦以及周边国际的经济文化交流有着重要意义。走廊穿越喀喇昆仑山系，在海拔4300-4750 m以上高寒山区，冻胀草丘，热融塘湖，石环、石链等冻土地貌广泛分布。正确的评价中巴走廊多年冻土分布问题是规划解决中巴经济走廊公路、铁路、通信和输油管线“四位一体”实际工程问题的基础，也对走廊内水资源利用、生态环境保护具有重要意义，尤其对边防建设有重要的战略意义。

传统的冻土调查通常采用野外工程勘查进行，但中巴走廊沿线地域广阔，实测资料较少，现有的多年冻土考察主要在中巴公路沿线，冻土资料基本一片空白。由于冻土与气候系统，特别是温度条件关系密切，国内外依据冻土与气温相互作用特性发展了一系列冻土分布模型。Smith等在分析气候-冻土关系基础上，考虑局部地形和土壤条件提出一种联系地温和气温的，半经验半物理的多年冻土顶板温度模型（TTOP）[5]。TTOP模型目前已应用于加拿大全国尺度和区域尺度，如Machensize河流域、布拉多高原和昂加瓦湾等的冻土分布模拟制图上，取得了良好的效果[6-8]。

随着遥感技术的快速发展，通过遥感技术可以大面积多时相地获得与多年冻土分布有关的地面信息进行空间分析，结合现有的冻土空间分布模型，可以实现宏观、动态、快速的大尺度上的冻土分布模拟。本研究以中巴经济走廊为研究区，包括中国新疆的喀什地区、克孜勒苏柯尔克孜自治州以及巴基斯坦全境（图1）。利用2013-2017年MODIS地表温度产品作为数据源代替实测数据输入到TTOP模型中，制备中巴走廊多年冻土分布数据集。

1 数据采集和处理方法

1.1 数据来源与预处理

中巴经济走廊多年冻土分布数据制备基于以下数据：2013-2017年 MODIS地表温度产品MOD11A1 LST与MYD11A1LST数据，2009年中国帕米尔高原第二次冰川编目数据集（v1.0），2003-2004年巴基斯坦冰川编目数据集，2008年世界土壤数据库（Harmonized World Soil Database version 1.2，即HWSDv1.2），如表1所示。

表1 研究采用数据列表

地表温度作为地表能量平衡中的主要参数以及气候系统的主要因子，能够表征地气间能量和水分交换的程度，是影响冻土发育、分布以及演化的关键因素之一，也是冻土建模的上边界条件。MOD11A1与MYD11A1均采用劈窗法从MOIDS数据中反演得到地表温度，为每日的L3级产品，空间分辨率1 km，均包括白天和夜间数据，来源于美国国家航空航天局（NASA）陆地过程分布式数据档案中心（LPDAAC）。数据以HDF-EOS格式存储，并采用正弦曲线投影坐标系。覆盖整个中巴走廊需要6景数据。

图1 研究区范围地理位置示意图

本研究获取 2013-2017年 MOD11A1白天地表温度数据、MOD11A1夜间地表温度数据、MYD11A1白天地表温度数据及 MYD11A1夜间地表温度数据等 4种数据产品共 21 900景。使用MRT工具对数据进行镶嵌和重采样，将所有数据坐标系转为WGS 84坐标系，再利用中巴走廊矢量边界和开源栅格处理库GDAL中的warp工具进行裁剪。

中国帕米尔高原第二次冰川编目数据集和巴基斯坦冰川编目均采用Landsat TM/ETM+数据经过校正后通过自动提取和专家干预修订获得，数据来源于寒区旱区科学数据中心（WestDC）[9-10]，将两种数据矢量合并后依照研究区范围进行裁剪。

TTOP模型需要不同类型土壤的冻融热导系数作为输入参数。土壤数据来源于2008年维也纳国际应用系统分析研究所（IIASA）和联合国粮农组织（FAO）所构建的HWSDv1.2。HWSDv1.2包括不同种类土壤的全球分布情况，土壤分类标准为FAO 90，数据分辨率1 km，数据投影WGS 84。

将冰川编目数据和土壤类型数据使用中巴经济走廊矢量边界裁剪，形成研究区内的冰川分布数据和土壤类型分布数据。

1.2 数据处理步骤

1.2.1 TTOP模型

TTOP模型其实质是采用活动层顶板温度来估计多年冻土的热状况。模型假定土壤基质匀质，在年尺度上，地-气热交换界面处于热平衡状态。TTOP模型定义如下：

式中，TTOP为多年冻土判别指标，TTOP<0表示存在多年冻土，TTOP>0表示季节性冻土或无冻土；P为周期，取365天；Kt和Kf分别为土壤融化和冻结时的导热系数；DDT为融化指数，指一年中连续高于0℃的持续时间和其数值乘积的总和，以℃∙day表示；DDF为冻结指数，指一年中低于0℃的持续时间和其数值乘积的总和，以℃∙day表示。nt和nf分别为融化和冻结因子，是由于地表植被和雪盖会影响地气热交换而引入融化和冻结因子修正模型结果。中巴走廊多年冻土发育地区处于高寒山区，也是极端干旱区，植被稀疏，本研究忽略积雪覆盖和植被影响，将nt和nf都赋值为 1。

1.2.2 地表温度数据插值

MODIS LST数据每一景在研究区内都有较多空缺值，需在时间序列上对空缺值进行插值。地表温度通常表现出明显的周期性，李述训等采用余弦函数对青藏高原时间序列地表温度进行模拟[11]，建立了地表温度与时间的函数关系，表达式如下：

其中，Ts为地表温度，t为时间，为时间序列内的平均温度；A为时间序列内的温度振幅；ω=2π/T中，T为时间序列周期，φ为相位角。

上述模拟方法虽能够反映地表温度周期变化情况，但参数计算不够准确。以平均温度和温度振幅作为余弦函数系数使函数拟合后整体偏差较大。本研究以T̅、A和φ为待定系数，利用最小二乘法对全年365景数据每一景数据对应像元进行拟合，利用拟合函数对地表温度数据的空缺值做插值处理。

1.2.3 导热系数

土壤融化和冻结时的导热系数是TTOP模型中与土壤性质有关的参数，其值取决于土壤颗粒级配、干容重、矿物成分和含水量等因素。由于缺少中巴走廊土壤性质的详细参数，难以计算每一网格上的土壤导热系数。王之夏等在研究青藏高原多年冻土分布时基于青藏高原第四纪地图，将土壤按地质成因分类，然后针对每一类采用Johansen方法计算其冻融状态下的导热系数[12-14]。由于此方法基于土壤本身性质，地区差异性相对较小，因此本研究认为可用在中巴经济走廊。HWSD数据包含土壤分类信息，可获取栅格上土壤分类名称，但采用的FAO土壤分类体系是以土壤诊断层和诊断特性为基础，经过对比筛选后得到基于HWSD数据的冻融导热系数，见表2。

表2 不同类型土壤的冻融导热系数表

1.2.4 冰川影响

本研究不考虑冰川下多年冻土情况。以地表温度和TTOP模型计算冻土分布会将部分冰川错分成冻土，需要根据研究区内冰川分布数据将其去除。

1.2.5 多年冻土分布模拟流程

首先将2013-2017每年MOD11A1与MYD11A1两种数据集的白天地表温度数据LST_Day_1km做算术平均，夜晚地表温度数据LST_Night_1km做算术平均；然后将每日白天地表温度平均值与夜晚地表温度平均值做融合，得到每日平均温度。融合方式为：白天均值与夜晚均值都存在则做算术平均得到当日平均温度值，有一个值空缺则当日平均温度值空缺；再做地表温度数据时间序列插值，将5a中相应日期的5个数据进行平均后计算DDT和DDF，制作研究区土壤融化导热系数图和土壤冻结导热系数图；最后将以上数据输入TTOP模型中，剔除冰川覆盖区域，通过二值分类得到中巴走廊多年冻土分布图。数据处理流程如图2所示：

图2 冻土分布数据处理流程图

2 数据样本描述

本数据集分为2个文件夹存储，分别为冻土分布栅格数据和冻土分布矢量数据。前者文件夹包括冻土分布栅格图CPEC_PermafrostDistribution，空间分辨率1 km，保存格式为tif。后者文件夹内为多年冻土分布矢量图CPEC_PermafrostDistribution，采用SHP格式保存。所有数据地理坐标系为WGS 84。数据结果展示如图3。

3 数据质量控制和评估

3.1 数据重建质量评估

以余弦函数和最小二乘法拟合的 MODIS地表温度数据可用决定系数来评价。决定系数表示一个随机变量与多个随机变量关系的数字特征，用来反映模型对样本数据的近似程度，常以R2表示。R2值域为（0,1），越接近1表明拟合的优度越高，其表达式如下：

图3 中巴经济走廊多年冻土分布图

中巴经济走廊边界区域内共有1 238 512像素点，将2013-2017年日平均地表温度做时间序列拟合，其决定系数R2最大值为0.95185，最小值0.0014。其空间分布情况与频率分布直方图如图4所示。由图可知，中巴走廊大部分区域整体拟合优度较高（R2≥0.6），拟合优度较差（R2＜0.4）区域则主要出现在帕米尔高原冰川和巴基斯坦沿海。而本研究不考虑冰川覆盖下冻土分布情况，同时巴基斯坦沿海气候与海拔不具备冻土形成条件，因此上述拟合优度较差区域对冻土模拟分布影响不大。所以，时间序列拟合的方式在中巴走廊冻土分布模拟研究中是可用的。

图4 决定系数R2空间分布图与频率分布直方图

3.2 多年冻土分布模拟结果评估

由于缺少冻土分布的实地勘测资料，因此以文献所提中巴公路沿线冻土数据作为验证。张学进等人提出中巴多年冻土分布于海拔4500 m以上的红其拉普路段K807+000至K811+343.165路段内，连续分布长度4.3 km[15]。朱颖彦等人提到中巴走廊多年冻土分布在红其拉普附近，海拔4300 m以上路段[16]。两篇文献提到多年冻土位置接近，可互为验证。本研究以遥感数据和模型模拟的方式得到多年冻土分布图，在中巴公路沿线所得结果如图5所示，与以上两篇文献吻合。

图5 红其拉普冻土分布图

4 数据价值

传统冻土图绘制主要以野外实测数据为基础，如钻孔地温数据、常规气象数据等。虽然结果精度较高，但其必须建立在时间序列足够长、观测数据足够可靠的基础上。中巴经济走廊多年冻土分布区域自然环境严苛，地缘政治环境复杂，在此区域进行冻土的调查和原位观测具有较多限制。

本文基于遥感数据MODIS地表温度数据产品，利用TTOP模型模拟中巴走廊冻土分布情况。针对数据与模型时间尺度不匹配的情况，采用拟合方式对数据进行时间序列重建。使用世界土壤数据集数据估算不同种类土壤冻融导热系数作为TTOP模型参数，最后获得走廊多年冻土分布图，并采用决定系数对数据重建过程中的拟合优度进行评价。野外实测是调查冻土分布的最佳手段，本文搜集了文献中关于多年冻土分布位置的描述，通过和模拟结果进行比较可知，基于MODIS地表温度数据和TTOP模型的多年冻土分布模拟结果与文献中的野外实测结果具有较高的一致性。

中巴经济走廊多年冻土分布数据是走廊内工程选址、设计和规划的重要基础数据，也可作为全球气候变化背景下，走廊多年冻土变化的本底数据。在多年冻土分布数据基础上，结合气候、水文、生态等数据综合分析，对区域内工程长期安全运营和生态环境可持续发展有着重要意义。

5 数据使用方法和建议

2017年中巴经济走廊冻土分布数据集保存为栅格TIF格式和矢量SHP格式。ArcGIS、QGIS、ENVI、ERDAS等常用的GIS与遥感软件可支持该数据的读取和操作。

致 谢

感谢LPDAAC提供的MODIS数据。感谢FAO提供的世界土壤数据库HWSD v1.2。感谢寒区旱区科学数据中心提供的中国帕米尔高原第二次冰川编目数据和巴基斯坦冰川编目数据。感谢中科院超级计算兰州分中心提供计算资源和存储资源。