2017年中巴經(jīng)濟(jì)走廊凍土分布數(shù)據(jù)集

艾鳴浩，張耀南＊，康建芳，馮克庭，田德宇

1. 中國(guó)科學(xué)院西北生態(tài)環(huán)境研究院，科學(xué)大數(shù)據(jù)中心，蘭州 730070

2. 國(guó)家特殊環(huán)境、特殊功能觀測(cè)研究臺(tái)站共享服務(wù)平臺(tái)，平臺(tái)服務(wù)中心，蘭州730070

數(shù)據(jù)庫(kù)（集）基本信息簡(jiǎn)介

引 言

凍土是北半球冰凍圈分布最廣泛的因子，占陸地面積的56%，其中多年凍土約占24%[1]。多年凍土具有熱穩(wěn)定性差、水熱活動(dòng)強(qiáng)烈、厚層地下冰和高含冰量?jī)鐾了急戎卮蟮忍攸c(diǎn)，在全球氣候變暖背景下，多年凍土對(duì)工程活動(dòng)的熱擾動(dòng)極為敏感。在多年凍土區(qū)工程建設(shè)中，施工期不可避免地會(huì)改變包括地表覆蓋、地形因子及淺層地表能量平衡等局地因素，加上開(kāi)通運(yùn)營(yíng)后動(dòng)載荷對(duì)凍融循環(huán)水熱環(huán)境的改變，直接導(dǎo)致凍土溫度升高和上限抬升，進(jìn)而增加了凍脹、融沉和凍拔等凍土災(zāi)害風(fēng)險(xiǎn)，對(duì)多年凍土區(qū)工程施工和安全運(yùn)營(yíng)構(gòu)成嚴(yán)重威脅，也對(duì)區(qū)域水文地質(zhì)條件和生態(tài)環(huán)境演化產(chǎn)生巨大影響[2-4]。

中巴經(jīng)濟(jì)走廊作為“一帶一路”倡議的重要組成部分，是中國(guó)與中亞、西亞、南亞地區(qū)的交通要道，對(duì)中國(guó)和巴基斯坦以及周邊國(guó)際的經(jīng)濟(jì)文化交流有著重要意義。走廊穿越喀喇昆侖山系，在海拔4300-4750 m以上高寒山區(qū)，凍脹草丘，熱融塘湖，石環(huán)、石鏈等凍土地貌廣泛分布。正確的評(píng)價(jià)中巴走廊多年凍土分布問(wèn)題是規(guī)劃解決中巴經(jīng)濟(jì)走廊公路、鐵路、通信和輸油管線“四位一體”實(shí)際工程問(wèn)題的基礎(chǔ)，也對(duì)走廊內(nèi)水資源利用、生態(tài)環(huán)境保護(hù)具有重要意義，尤其對(duì)邊防建設(shè)有重要的戰(zhàn)略意義。

傳統(tǒng)的凍土調(diào)查通常采用野外工程勘查進(jìn)行，但中巴走廊沿線地域廣闊，實(shí)測(cè)資料較少，現(xiàn)有的多年凍土考察主要在中巴公路沿線，凍土資料基本一片空白。由于凍土與氣候系統(tǒng)，特別是溫度條件關(guān)系密切，國(guó)內(nèi)外依據(jù)凍土與氣溫相互作用特性發(fā)展了一系列凍土分布模型。Smith等在分析氣候-凍土關(guān)系基礎(chǔ)上，考慮局部地形和土壤條件提出一種聯(lián)系地溫和氣溫的，半經(jīng)驗(yàn)半物理的多年凍土頂板溫度模型（TTOP）[5]。TTOP模型目前已應(yīng)用于加拿大全國(guó)尺度和區(qū)域尺度，如Machensize河流域、布拉多高原和昂加瓦灣等的凍土分布模擬制圖上，取得了良好的效果[6-8]。

隨著遙感技術(shù)的快速發(fā)展，通過(guò)遙感技術(shù)可以大面積多時(shí)相地獲得與多年凍土分布有關(guān)的地面信息進(jìn)行空間分析，結(jié)合現(xiàn)有的凍土空間分布模型，可以實(shí)現(xiàn)宏觀、動(dòng)態(tài)、快速的大尺度上的凍土分布模擬。本研究以中巴經(jīng)濟(jì)走廊為研究區(qū)，包括中國(guó)新疆的喀什地區(qū)、克孜勒蘇柯?tīng)柨俗巫灾沃菀约鞍突固谷常▓D1）。利用2013-2017年MODIS地表溫度產(chǎn)品作為數(shù)據(jù)源代替實(shí)測(cè)數(shù)據(jù)輸入到TTOP模型中，制備中巴走廊多年凍土分布數(shù)據(jù)集。

1 數(shù)據(jù)采集和處理方法

1.1 數(shù)據(jù)來(lái)源與預(yù)處理

中巴經(jīng)濟(jì)走廊多年凍土分布數(shù)據(jù)制備基于以下數(shù)據(jù)：2013-2017年 MODIS地表溫度產(chǎn)品MOD11A1 LST與MYD11A1LST數(shù)據(jù)，2009年中國(guó)帕米爾高原第二次冰川編目數(shù)據(jù)集（v1.0），2003-2004年巴基斯坦冰川編目數(shù)據(jù)集，2008年世界土壤數(shù)據(jù)庫(kù)（Harmonized World Soil Database version 1.2，即HWSDv1.2），如表1所示。

表1 研究采用數(shù)據(jù)列表

地表溫度作為地表能量平衡中的主要參數(shù)以及氣候系統(tǒng)的主要因子，能夠表征地氣間能量和水分交換的程度，是影響凍土發(fā)育、分布以及演化的關(guān)鍵因素之一，也是凍土建模的上邊界條件。MOD11A1與MYD11A1均采用劈窗法從MOIDS數(shù)據(jù)中反演得到地表溫度，為每日的L3級(jí)產(chǎn)品，空間分辨率1 km，均包括白天和夜間數(shù)據(jù)，來(lái)源于美國(guó)國(guó)家航空航天局（NASA）陸地過(guò)程分布式數(shù)據(jù)檔案中心（LPDAAC）。數(shù)據(jù)以HDF-EOS格式存儲(chǔ)，并采用正弦曲線投影坐標(biāo)系。覆蓋整個(gè)中巴走廊需要6景數(shù)據(jù)。

圖1 研究區(qū)范圍地理位置示意圖

本研究獲取 2013-2017年 MOD11A1白天地表溫度數(shù)據(jù)、MOD11A1夜間地表溫度數(shù)據(jù)、MYD11A1白天地表溫度數(shù)據(jù)及 MYD11A1夜間地表溫度數(shù)據(jù)等 4種數(shù)據(jù)產(chǎn)品共 21 900景。使用MRT工具對(duì)數(shù)據(jù)進(jìn)行鑲嵌和重采樣，將所有數(shù)據(jù)坐標(biāo)系轉(zhuǎn)為WGS 84坐標(biāo)系，再利用中巴走廊矢量邊界和開(kāi)源柵格處理庫(kù)GDAL中的warp工具進(jìn)行裁剪。

中國(guó)帕米爾高原第二次冰川編目數(shù)據(jù)集和巴基斯坦冰川編目均采用Landsat TM/ETM+數(shù)據(jù)經(jīng)過(guò)校正后通過(guò)自動(dòng)提取和專(zhuān)家干預(yù)修訂獲得，數(shù)據(jù)來(lái)源于寒區(qū)旱區(qū)科學(xué)數(shù)據(jù)中心（WestDC）[9-10]，將兩種數(shù)據(jù)矢量合并后依照研究區(qū)范圍進(jìn)行裁剪。

TTOP模型需要不同類(lèi)型土壤的凍融熱導(dǎo)系數(shù)作為輸入?yún)?shù)。土壤數(shù)據(jù)來(lái)源于2008年維也納國(guó)際應(yīng)用系統(tǒng)分析研究所（IIASA）和聯(lián)合國(guó)糧農(nóng)組織（FAO）所構(gòu)建的HWSDv1.2。HWSDv1.2包括不同種類(lèi)土壤的全球分布情況，土壤分類(lèi)標(biāo)準(zhǔn)為FAO 90，數(shù)據(jù)分辨率1 km，數(shù)據(jù)投影WGS 84。

將冰川編目數(shù)據(jù)和土壤類(lèi)型數(shù)據(jù)使用中巴經(jīng)濟(jì)走廊矢量邊界裁剪，形成研究區(qū)內(nèi)的冰川分布數(shù)據(jù)和土壤類(lèi)型分布數(shù)據(jù)。

1.2 數(shù)據(jù)處理步驟

1.2.1 TTOP模型

TTOP模型其實(shí)質(zhì)是采用活動(dòng)層頂板溫度來(lái)估計(jì)多年凍土的熱狀況。模型假定土壤基質(zhì)勻質(zhì)，在年尺度上，地-氣熱交換界面處于熱平衡狀態(tài)。TTOP模型定義如下：

式中，TTOP為多年凍土判別指標(biāo)，TTOP<0表示存在多年凍土，TTOP>0表示季節(jié)性凍土或無(wú)凍土；P為周期，取365天；Kt和Kf分別為土壤融化和凍結(jié)時(shí)的導(dǎo)熱系數(shù)；DDT為融化指數(shù)，指一年中連續(xù)高于0℃的持續(xù)時(shí)間和其數(shù)值乘積的總和，以℃?day表示；DDF為凍結(jié)指數(shù)，指一年中低于0℃的持續(xù)時(shí)間和其數(shù)值乘積的總和，以℃?day表示。nt和nf分別為融化和凍結(jié)因子，是由于地表植被和雪蓋會(huì)影響地氣熱交換而引入融化和凍結(jié)因子修正模型結(jié)果。中巴走廊多年凍土發(fā)育地區(qū)處于高寒山區(qū)，也是極端干旱區(qū)，植被稀疏，本研究忽略積雪覆蓋和植被影響，將nt和nf都賦值為 1。

1.2.2 地表溫度數(shù)據(jù)插值

MODIS LST數(shù)據(jù)每一景在研究區(qū)內(nèi)都有較多空缺值，需在時(shí)間序列上對(duì)空缺值進(jìn)行插值。地表溫度通常表現(xiàn)出明顯的周期性，李述訓(xùn)等采用余弦函數(shù)對(duì)青藏高原時(shí)間序列地表溫度進(jìn)行模擬[11]，建立了地表溫度與時(shí)間的函數(shù)關(guān)系，表達(dá)式如下：

其中，Ts為地表溫度，t為時(shí)間，為時(shí)間序列內(nèi)的平均溫度；A為時(shí)間序列內(nèi)的溫度振幅；ω=2π/T中，T為時(shí)間序列周期，φ為相位角。

上述模擬方法雖能夠反映地表溫度周期變化情況，但參數(shù)計(jì)算不夠準(zhǔn)確。以平均溫度和溫度振幅作為余弦函數(shù)系數(shù)使函數(shù)擬合后整體偏差較大。本研究以T?、A和φ為待定系數(shù)，利用最小二乘法對(duì)全年365景數(shù)據(jù)每一景數(shù)據(jù)對(duì)應(yīng)像元進(jìn)行擬合，利用擬合函數(shù)對(duì)地表溫度數(shù)據(jù)的空缺值做插值處理。

1.2.3 導(dǎo)熱系數(shù)

土壤融化和凍結(jié)時(shí)的導(dǎo)熱系數(shù)是TTOP模型中與土壤性質(zhì)有關(guān)的參數(shù)，其值取決于土壤顆粒級(jí)配、干容重、礦物成分和含水量等因素。由于缺少中巴走廊土壤性質(zhì)的詳細(xì)參數(shù)，難以計(jì)算每一網(wǎng)格上的土壤導(dǎo)熱系數(shù)。王之夏等在研究青藏高原多年凍土分布時(shí)基于青藏高原第四紀(jì)地圖，將土壤按地質(zhì)成因分類(lèi)，然后針對(duì)每一類(lèi)采用Johansen方法計(jì)算其凍融狀態(tài)下的導(dǎo)熱系數(shù)[12-14]。由于此方法基于土壤本身性質(zhì)，地區(qū)差異性相對(duì)較小，因此本研究認(rèn)為可用在中巴經(jīng)濟(jì)走廊。HWSD數(shù)據(jù)包含土壤分類(lèi)信息，可獲取柵格上土壤分類(lèi)名稱(chēng)，但采用的FAO土壤分類(lèi)體系是以土壤診斷層和診斷特性為基礎(chǔ)，經(jīng)過(guò)對(duì)比篩選后得到基于HWSD數(shù)據(jù)的凍融導(dǎo)熱系數(shù)，見(jiàn)表2。

表2 不同類(lèi)型土壤的凍融導(dǎo)熱系數(shù)表

1.2.4 冰川影響

本研究不考慮冰川下多年凍土情況。以地表溫度和TTOP模型計(jì)算凍土分布會(huì)將部分冰川錯(cuò)分成凍土，需要根據(jù)研究區(qū)內(nèi)冰川分布數(shù)據(jù)將其去除。

1.2.5 多年凍土分布模擬流程

首先將2013-2017每年MOD11A1與MYD11A1兩種數(shù)據(jù)集的白天地表溫度數(shù)據(jù)LST_Day_1km做算術(shù)平均，夜晚地表溫度數(shù)據(jù)LST_Night_1km做算術(shù)平均；然后將每日白天地表溫度平均值與夜晚地表溫度平均值做融合，得到每日平均溫度。融合方式為：白天均值與夜晚均值都存在則做算術(shù)平均得到當(dāng)日平均溫度值，有一個(gè)值空缺則當(dāng)日平均溫度值空缺；再做地表溫度數(shù)據(jù)時(shí)間序列插值，將5a中相應(yīng)日期的5個(gè)數(shù)據(jù)進(jìn)行平均后計(jì)算DDT和DDF，制作研究區(qū)土壤融化導(dǎo)熱系數(shù)圖和土壤凍結(jié)導(dǎo)熱系數(shù)圖；最后將以上數(shù)據(jù)輸入TTOP模型中，剔除冰川覆蓋區(qū)域，通過(guò)二值分類(lèi)得到中巴走廊多年凍土分布圖。數(shù)據(jù)處理流程如圖2所示：

圖2 凍土分布數(shù)據(jù)處理流程圖

2 數(shù)據(jù)樣本描述

本數(shù)據(jù)集分為2個(gè)文件夾存儲(chǔ)，分別為凍土分布柵格數(shù)據(jù)和凍土分布矢量數(shù)據(jù)。前者文件夾包括凍土分布柵格圖CPEC_PermafrostDistribution，空間分辨率1 km，保存格式為tif。后者文件夾內(nèi)為多年凍土分布矢量圖CPEC_PermafrostDistribution，采用SHP格式保存。所有數(shù)據(jù)地理坐標(biāo)系為WGS 84。數(shù)據(jù)結(jié)果展示如圖3。

3 數(shù)據(jù)質(zhì)量控制和評(píng)估

3.1 數(shù)據(jù)重建質(zhì)量評(píng)估

以余弦函數(shù)和最小二乘法擬合的 MODIS地表溫度數(shù)據(jù)可用決定系數(shù)來(lái)評(píng)價(jià)。決定系數(shù)表示一個(gè)隨機(jī)變量與多個(gè)隨機(jī)變量關(guān)系的數(shù)字特征，用來(lái)反映模型對(duì)樣本數(shù)據(jù)的近似程度，常以R2表示。R2值域?yàn)椋?,1），越接近1表明擬合的優(yōu)度越高，其表達(dá)式如下：

圖3 中巴經(jīng)濟(jì)走廊多年凍土分布圖

中巴經(jīng)濟(jì)走廊邊界區(qū)域內(nèi)共有1 238 512像素點(diǎn)，將2013-2017年日平均地表溫度做時(shí)間序列擬合，其決定系數(shù)R2最大值為0.95185，最小值0.0014。其空間分布情況與頻率分布直方圖如圖4所示。由圖可知，中巴走廊大部分區(qū)域整體擬合優(yōu)度較高（R2≥0.6），擬合優(yōu)度較差（R2＜0.4）區(qū)域則主要出現(xiàn)在帕米爾高原冰川和巴基斯坦沿海。而本研究不考慮冰川覆蓋下凍土分布情況，同時(shí)巴基斯坦沿海氣候與海拔不具備凍土形成條件，因此上述擬合優(yōu)度較差區(qū)域?qū)鐾聊M分布影響不大。所以，時(shí)間序列擬合的方式在中巴走廊凍土分布模擬研究中是可用的。

圖4 決定系數(shù)R2空間分布圖與頻率分布直方圖

3.2 多年凍土分布模擬結(jié)果評(píng)估

由于缺少凍土分布的實(shí)地勘測(cè)資料，因此以文獻(xiàn)所提中巴公路沿線凍土數(shù)據(jù)作為驗(yàn)證。張學(xué)進(jìn)等人提出中巴多年凍土分布于海拔4500 m以上的紅其拉普路段K807+000至K811+343.165路段內(nèi)，連續(xù)分布長(zhǎng)度4.3 km[15]。朱穎彥等人提到中巴走廊多年凍土分布在紅其拉普附近，海拔4300 m以上路段[16]。兩篇文獻(xiàn)提到多年凍土位置接近，可互為驗(yàn)證。本研究以遙感數(shù)據(jù)和模型模擬的方式得到多年凍土分布圖，在中巴公路沿線所得結(jié)果如圖5所示，與以上兩篇文獻(xiàn)吻合。

圖5 紅其拉普凍土分布圖

4 數(shù)據(jù)價(jià)值

傳統(tǒng)凍土圖繪制主要以野外實(shí)測(cè)數(shù)據(jù)為基礎(chǔ)，如鉆孔地溫?cái)?shù)據(jù)、常規(guī)氣象數(shù)據(jù)等。雖然結(jié)果精度較高，但其必須建立在時(shí)間序列足夠長(zhǎng)、觀測(cè)數(shù)據(jù)足夠可靠的基礎(chǔ)上。中巴經(jīng)濟(jì)走廊多年凍土分布區(qū)域自然環(huán)境嚴(yán)苛，地緣政治環(huán)境復(fù)雜，在此區(qū)域進(jìn)行凍土的調(diào)查和原位觀測(cè)具有較多限制。

本文基于遙感數(shù)據(jù)MODIS地表溫度數(shù)據(jù)產(chǎn)品，利用TTOP模型模擬中巴走廊凍土分布情況。針對(duì)數(shù)據(jù)與模型時(shí)間尺度不匹配的情況，采用擬合方式對(duì)數(shù)據(jù)進(jìn)行時(shí)間序列重建。使用世界土壤數(shù)據(jù)集數(shù)據(jù)估算不同種類(lèi)土壤凍融導(dǎo)熱系數(shù)作為T(mén)TOP模型參數(shù)，最后獲得走廊多年凍土分布圖，并采用決定系數(shù)對(duì)數(shù)據(jù)重建過(guò)程中的擬合優(yōu)度進(jìn)行評(píng)價(jià)。野外實(shí)測(cè)是調(diào)查凍土分布的最佳手段，本文搜集了文獻(xiàn)中關(guān)于多年凍土分布位置的描述，通過(guò)和模擬結(jié)果進(jìn)行比較可知，基于MODIS地表溫度數(shù)據(jù)和TTOP模型的多年凍土分布模擬結(jié)果與文獻(xiàn)中的野外實(shí)測(cè)結(jié)果具有較高的一致性。

中巴經(jīng)濟(jì)走廊多年凍土分布數(shù)據(jù)是走廊內(nèi)工程選址、設(shè)計(jì)和規(guī)劃的重要基礎(chǔ)數(shù)據(jù)，也可作為全球氣候變化背景下，走廊多年凍土變化的本底數(shù)據(jù)。在多年凍土分布數(shù)據(jù)基礎(chǔ)上，結(jié)合氣候、水文、生態(tài)等數(shù)據(jù)綜合分析，對(duì)區(qū)域內(nèi)工程長(zhǎng)期安全運(yùn)營(yíng)和生態(tài)環(huán)境可持續(xù)發(fā)展有著重要意義。

5 數(shù)據(jù)使用方法和建議

2017年中巴經(jīng)濟(jì)走廊凍土分布數(shù)據(jù)集保存為柵格TIF格式和矢量SHP格式。ArcGIS、QGIS、ENVI、ERDAS等常用的GIS與遙感軟件可支持該數(shù)據(jù)的讀取和操作。

致 謝

感謝LPDAAC提供的MODIS數(shù)據(jù)。感謝FAO提供的世界土壤數(shù)據(jù)庫(kù)HWSD v1.2。感謝寒區(qū)旱區(qū)科學(xué)數(shù)據(jù)中心提供的中國(guó)帕米爾高原第二次冰川編目數(shù)據(jù)和巴基斯坦冰川編目數(shù)據(jù)。感謝中科院超級(jí)計(jì)算蘭州分中心提供計(jì)算資源和存儲(chǔ)資源。