黃博文，高瑞翔，陳一仰，范增輝，劉修國

(中國地質(zhì)大學(xué)(武漢)地理與信息工程學(xué)院，湖北 武漢430074)

凍土是指土壤溫度在零攝氏度及以下，并含有冰的各種巖石和土壤，一般可分為短時凍土(數(shù)小時/數(shù)日以至半月)、季節(jié)凍土(半月至數(shù)月)和多年凍土(又稱永久凍土，指的是持續(xù)2年或2年以上的 凍結(jié)不融的土層)。在全球變暖的影響下，凍土已經(jīng)成為全球氣候變化的重要觀測對象之一，通過對凍土的觀測可以一定程度上反映全球變暖的嚴(yán)重程度[1]。

青藏高原作為世界上中低緯度海拔最高、面積最大的多年凍土區(qū)域，其凍土變化會對中國東部乃 至東亞氣候的形成、變化和發(fā)展造成重大的影響。由于青藏高原海拔高、地域廣、環(huán)境惡劣，通過地面調(diào)查研究凍土環(huán)境的工作方法難以大范圍地應(yīng)用。許多基于溫度數(shù)據(jù)的凍土分布模型往往需要依靠地表溫度觀測站點獲取地表溫度數(shù)據(jù)，具有空間局限性明顯、監(jiān)測密度低、人力物力成本較高等特點，也難以在區(qū)域尺度上應(yīng)用。遙感(RS)技術(shù)能夠在一定的條件下進(jìn)行大范圍的數(shù)據(jù)采集，較好地解決了由于環(huán)境險惡造成的交通不便、測量困難等難題，彌補了傳統(tǒng)研究技術(shù)上的不足。本文針對青藏高原凍土的特點，提出一種基于多源RS數(shù)據(jù)融合技術(shù)反演地表溫度，并結(jié)合相關(guān)的算法，實現(xiàn)凍土空間分布提取和分析的方法。該方法可以消除極端天氣對地表溫度反演結(jié)果的不利影響，實現(xiàn)凍土環(huán)境的精細(xì)化監(jiān)測。

1 研究區(qū)概況與RS數(shù)據(jù)源

1.1 研究區(qū)概況

本文研究區(qū)為青藏高原地區(qū)。青藏高原位于北緯25°～40°E、東經(jīng)74°～104°E之間，面積約250萬km2，平均海拔為4 500 m。主要地形地貌為山地，其中包括東西與南北走向的兩種山脈，其間還分布大量的湖泊與河流。該地區(qū)氣溫常年較低，部分地區(qū)年積溫低于0℃，為凍土形成創(chuàng)造了先決的條件。 青藏高原多年凍土區(qū)面積約為157萬km2,主要以亞穩(wěn)定型(30.4%)、過渡型(22.1%)和不穩(wěn)定型(22.6%)多年凍土為主，極穩(wěn)定型和穩(wěn)定型多年凍土的占比較少。此外，選擇黑山實驗區(qū)為重點研究區(qū)。它位于凍土北界，庫塞湖東北方向，在柴達(dá)木盆地東北邊緣，深藏于祁連山腹地。該實驗區(qū)毗鄰黑海，地形多樣,地表覆蓋類型眾多,包括植被、冰川、山地、湖泊等，有利于檢驗RS監(jiān)測中多源數(shù)據(jù)快速融合并提取凍土的方法。

1.2 RS數(shù)據(jù)源的選擇及其特征

目前國內(nèi)外提出了許多基于溫度數(shù)據(jù)的凍土分布模型如凍結(jié)指數(shù)模型[2]，還有利用SAR數(shù)據(jù)進(jìn)行凍土地表覆蓋提取[3]等。隨著遙感科學(xué)的發(fā)展，通過遙感影像反演地表溫度的技術(shù)已經(jīng)較為成熟，中分辨率成像光譜儀(MODIS)溫度產(chǎn)品數(shù)據(jù)已在凍土監(jiān)測研究中得到了應(yīng)用[4]。趙全寧等[5]通過線性回歸分析指出，在未來氣候變暖背景下，在玉樹地區(qū)應(yīng)用溫度因子估算平均最大凍土深度的變化具有較高的可信度；秦艷慧等[6]對ERA-Interim進(jìn)行了重新擬合改良,并使用其通過擬合得到了我國青藏高原地區(qū)的凍土分布圖，與實際情況對比其擬合效果較好。此外，近年來的相關(guān)研究充分說明了使用遙感相關(guān)的溫度產(chǎn)品反演青藏高原地區(qū)的凍土分布是可行的。

本文選擇了MODIS的8天合成地表溫度數(shù)據(jù)(MOD11A2)作為大尺度溫度數(shù)據(jù)，該數(shù)據(jù)時間分辨率高，涵蓋白天和夜晚地表溫度信息，采用的是8天內(nèi)晴好天氣狀況下的地表溫度數(shù)據(jù)，并根據(jù)地表覆蓋類型數(shù)據(jù)進(jìn)行平均計算而得到，能夠較好地反映地表在一段時間內(nèi)的溫度分布情況,其空間分辨率為1 km。依據(jù)該數(shù)據(jù)可以計算出相應(yīng)的年平均地溫和年積溫，用于凍土分布信息的自動化提取。小尺度影像選擇Landsat8 影像的2級數(shù)據(jù)，該等級數(shù)據(jù)經(jīng)過輻射校正和幾何校正，空間分辨率為30 m，重返周期約為16 d,圖像的精細(xì)度較高，可以利用紅外波段數(shù)據(jù)反演得到地表溫度結(jié)果。此外，水汽數(shù)據(jù)來源于MOD05-L2系列產(chǎn)品,其反映的是晴空區(qū)域陸地上空的云層上方的大氣柱水蒸氣量，通過該數(shù)據(jù)以及Landsat8的影像數(shù)據(jù)即可根據(jù)算法得到地表溫度。

地表覆蓋類型數(shù)據(jù)采用:①清華大學(xué)制作的全球地表覆蓋類型數(shù)據(jù)，該數(shù)據(jù)由清華大學(xué)宮鵬教授團(tuán)隊制作，包括2017年全球0.000 25度(對應(yīng)約30 m)土地覆蓋空間分布，包含農(nóng)田、森林、濕地、水體 等10類地表覆蓋類型；②美國NASA制作的全球1 km地表覆蓋類型分類。這兩種產(chǎn)品數(shù)據(jù)精度較高、分辨率尚可，經(jīng)過裁剪處理后可以作為中間變量直接使用。

2 地表溫度數(shù)據(jù)時空融合方法

2.1 數(shù)據(jù)預(yù)處理

2.1.1 MODIS數(shù)據(jù)預(yù)處理

青藏高原區(qū)域涉及6景MODIS數(shù)據(jù)，經(jīng)拼接裁剪后方可得到研究區(qū)的完整影像。全年的MOD11A2數(shù)據(jù)共計276景影像。本研究的時間跨度為2000年到2018年,計19年，共需要處理5 244景影像。使用IDL語言調(diào)用MCTK插件對MODIS數(shù)據(jù)進(jìn)行重投影、裁剪、鑲嵌等批量化預(yù)處理，較快地得到了MOD11溫度產(chǎn)品數(shù)據(jù)集。

由于大氣影響和云霧遮擋，MOD11A2數(shù)據(jù)中存在著異常值和無效值等隨機(jī)誤差，可依據(jù)拉依達(dá)準(zhǔn)則[7]進(jìn)行誤差消除。另外，由于地表溫度受太陽直射的影響大,僅使用日溫度數(shù)據(jù)將會導(dǎo)致反演的年平均地溫過高，因此需綜合使用MOD11A2的日夜溫度數(shù)據(jù)，采用日夜溫度平均的方式反演可信的地表溫度數(shù)據(jù)[4]。

2.1.2 Landsat影像數(shù)據(jù)預(yù)處理

Landsat8數(shù)據(jù)的預(yù)處理包括：裁剪、拼接、輻射定標(biāo)和去云處理。使用ENVI軟件中的內(nèi)嵌工具即可完成Landsat8數(shù)據(jù)的預(yù)處理。對于有云的部分，利用Landsat8數(shù)據(jù)提供的云層數(shù)據(jù)，對影像進(jìn)行掩膜處理，可以獲取去云后的 Landsat8數(shù)據(jù)。同時，可根據(jù)該數(shù)據(jù)提供的反射率去除高海拔的冰川，可避免冰川被誤分類而影響結(jié)果精度。

2.1.3 地表覆蓋類型數(shù)據(jù)預(yù)處理

對于下載的地表覆蓋類型數(shù)據(jù)文件，進(jìn)行拼接、裁剪處理，可獲得研究區(qū)內(nèi)的地物分類。該數(shù)據(jù)的主要作用是確定不同地表覆蓋類型下的土壤導(dǎo)熱系數(shù)，以提高地表溫度反演與凍土類型確定的計算精度。

2.2 Landsat與MODIS地表溫度數(shù)據(jù)的融合方法

2.2.1 Landsat地表溫度反演

傳統(tǒng)的單通道算法[8]是通過估計大氣參數(shù)修正單一通道中的大氣效應(yīng)來反演地表溫度。Practical Single-Channel(PSC)算法[9]在單通道算法的基礎(chǔ)上，針對不同傳感器計算固定參數(shù)，降低了單窗算法中由于普朗克函數(shù)線性化和大氣校正帶來的誤差。該算法分兩步進(jìn)行地表溫度TS估計：首先通過消除地表發(fā)射率和大氣效應(yīng)引起的誤差，計算地表黑體輻射亮溫B(TS)；然后根據(jù)普朗克函數(shù)的反函數(shù)，計算地表溫度TS。具體實現(xiàn)過程為通過公式(1)、(2)計算地表比輻射率，再通過公式(3)計算地表黑體輻射亮溫，將計算得到的地表黑體輻射亮溫代入公式(4)即可計算地表溫度結(jié)果。詳見公式(1)～(4)：

(1)

(2)

(3)

(4)

式中：PV表示植被覆蓋度，其中NDVI表示歸一化差異植被指數(shù)；NDVIS表示裸地的歸一化差異植被指數(shù)(取值為0.2)；NDVIV表示純植被的歸一化差異植被指數(shù)(取值為0.55);w表示水汽數(shù)據(jù)，來源于MOD05數(shù)據(jù);ε表示地表比輻射率；εS表示裸地的比輻射率(取值為0.966 8)；εV表示純植被的比輻射率(取值為0.986 3)；Lsen表示星上輻射亮度，為采用Landsat8的熱紅外波段輻射定標(biāo)之后的結(jié)果;參數(shù)a0～a7均為與傳感器有關(guān)的固定參數(shù)，本文使用Landsat配套參數(shù)，具體見表1；λ表示中心波長(nm)；C1= 1.191 04×108W·μm4/(m2·sr),C2=1.438 77×104μm·K。

表1 Landsat配套參數(shù)表Table 1 Parameter-value of Landsat

Landsat遙感影像在研究區(qū)的成像時間為中午時分，依照此方法反演出的地表溫度要比MOD11A2的日夜平均溫度高，使得最終的年平均地溫的計算結(jié)果較使用MOD11A2數(shù)據(jù)高出許多，是一個極大的不可忽略的誤差因素。不同RS數(shù)據(jù)源獲取的對應(yīng)時間段內(nèi)的地面溫度相互關(guān)聯(lián)，可以通過擬合分析技術(shù)研究其相關(guān)性[10]。本文采用了多項式擬合技術(shù)，將Landsat遙感影像數(shù)據(jù)反演出的地表溫度與對應(yīng)的MOD11A2地表溫度數(shù)據(jù)進(jìn)行擬合,MOD11A2數(shù)據(jù)的時相與Landsat遙感影像相一致，使得擬合結(jié)果更為可靠。采用三次多項式擬合，擬合結(jié)果即為地面當(dāng)日平均溫度，將擬合后的夏季與冬季溫度作為年最高與最低溫度計算年平均地溫。

2.2.2 凍融指數(shù)計算

凍融指數(shù)是凍土判識與類型區(qū)劃的依據(jù)，本文采用的計算方法與文獻(xiàn)[11]相同，其關(guān)鍵參數(shù)是確定累計正溫DDT和累計負(fù)溫DDF。基于地溫全年變化符合余弦曲線規(guī)律的假定，DDT和DDF的計算公式如下：

(6)

(7)

(8)

(9)

(10)

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

根據(jù)上述公式，只需使用全年極值溫度進(jìn)行計算即可，具體計算示意圖見圖1。此種處理方式的優(yōu)點在于數(shù)據(jù)處理量小，一年僅需要兩個極值月份的影像即可以計算凍融指數(shù)。根據(jù)凍土北界黑山實驗區(qū)內(nèi)五道梁站點的地表溫度觀測數(shù)據(jù)[12]，區(qū)內(nèi)最低溫度所在月份為1月或12月，最高溫度所在月份為7月或8月。在需要精細(xì)化凍土分布的區(qū)域，本文利用該計算方法對擬合處理后的Landsat溫度數(shù)據(jù)進(jìn)行凍融指數(shù)的估算。

圖1 地表溫度擬合曲線示意圖Fig.1 Sketch of surface temperature fitting curve 圖中：A為擬合曲線的振幅;Tp為擬合曲線的平均值所在處，可以參考圖中所標(biāo)注信息；DDT為累計正溫(℃)；DDF為累計負(fù)溫(℃)。

由于全球氣候變化的影響，氣候波動和極端氣候出現(xiàn)的次數(shù)有所增加，使得年內(nèi)溫度變化不再符合簡單的簡諧波動假設(shè)，會造成積溫計算的不準(zhǔn)確。 為此，本文提出一種針對MODIS全年溫度產(chǎn)品數(shù)據(jù)集的凍融指數(shù)計算方法，以0℃為界限，將一年的溫度劃分為正溫和負(fù)溫，并分別累加計算正積溫和負(fù)積溫,從而得到準(zhǔn)確的凍結(jié)指數(shù)和融化指數(shù)[13]。該方法從根本上消除了氣候波動對年內(nèi)積溫計算準(zhǔn)確性的影響，實驗表明該方法計算得到的年平均地溫的精度更高，有利于后續(xù)的凍土分類與制圖。

2.3 基于溫度的凍土類型識別與分類

以多源RS數(shù)據(jù)為基礎(chǔ)，本文利用TTOP(Temperature at the Top Of Permafrost)模型反演活動層頂板的溫度(TTOP)，其計算公式如下[11]：

(11)

式中：rk為凍融指數(shù)；p為土壤導(dǎo)熱系數(shù)[W/(m·K)]。

根據(jù)公式(11)，通過輸入土壤導(dǎo)熱系數(shù)和凍融指數(shù)，即可反演得到活動層頂板的溫度。利用反演出的活動層頂板溫度，根據(jù)程國棟院士提出的熱力學(xué)分布方法[14]，可以利用活動層頂板溫度實現(xiàn)凍土分類與制圖。

本文選取MODIS溫度產(chǎn)品2001年到2018年共18 年的數(shù)據(jù)，根據(jù)TTOP模型反演得到每年活動層頂板的平均溫度，對所有數(shù)據(jù)在時間域上進(jìn)行濾波處理,以降低異常年份對凍土反演結(jié)果的影響。將數(shù)據(jù)分為兩期，在每期結(jié)果內(nèi)對反演得到的活動層頂板溫度進(jìn)行平均，據(jù)9年平均結(jié)果得到最終的TTOP作為分類依據(jù)，見表2。

表2 凍土類型劃分標(biāo)準(zhǔn)Table 2 Permafrost classification standards

3 地表溫度數(shù)據(jù)融合和凍土類型RS識別

3.1 Landsat與MODIS地溫數(shù)據(jù)融合結(jié)果

利用上述Landsat與MODIS數(shù)據(jù)的融合方法，可得到2000—2018年青藏高原Landsat8地溫反演結(jié)果。本文任選2013年和2018年青藏高原Landsat8地表溫度的反演結(jié)果進(jìn)行分析，見圖2和圖3。

圖2 2013年青藏高原Landsat8地表溫度的反演結(jié)果Fig.2 Inversion result of surface temperature of Qinghai- Tibet Plateau in 2013 based on Landsat8

圖3 2018年青藏高原Landsat8地表溫度的反演結(jié)果Fig.3 Inversion result of surface temperature of Qinghai- Tibet Plateau in 2018 based on Landsat8

3.2 基于地表溫度的凍土類型RS識別結(jié)果

多年年平均活動層頂板溫度可由活動層頂板溫度TTOP數(shù)據(jù)經(jīng)平均化處理后得到，進(jìn)而反演得到青藏高原的凍土分布圖，見圖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 與歷史數(shù)據(jù)的一致性分析

將2001—2009年青藏高原凍土分布圖(見圖4)和2010—2018年青藏高原凍土分布圖(見圖5)分別與1∶300萬青藏高原凍土分布圖(見圖6)矢量數(shù)據(jù)和青藏高原歷史凍土分布圖(見圖7)進(jìn)行宏觀上的疊加比對，其結(jié)果見圖8和圖9。結(jié)果發(fā)現(xiàn)：1∶300萬青藏高原凍土的分布范圍與圖4和圖5基本一致，說明在宏觀上采用本文的凍土環(huán)境RS信息提取的技術(shù)方法是可信的[15]；在多年凍土的集中分布區(qū)域(如圖6中箭頭所示)，青藏高原歷史凍土分布(見圖7)與圖4、圖5基本吻合;2001—2009年青藏高原凍土分布與歷史凍土分布的相似性達(dá)90.42%(見圖8)，2010—2018年青藏高原凍土分布與歷史凍土分布的相似性達(dá)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的年平均地溫擬合的準(zhǔn)確度

以青藏高原2013年和2018年度為例，將Landsat反演得到的年平均地溫數(shù)據(jù)與MODIS標(biāo)準(zhǔn)溫度進(jìn)行了對比，其差值結(jié)果見圖12和圖13。

圖12 利用Landsat反演的2013年青藏高原年平均地 溫與MODIS標(biāo)準(zhǔn)溫度的差值分布圖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標(biāo)準(zhǔn)溫度的差值分布圖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標(biāo)準(zhǔn)溫度的平均誤差為0.419，方差為3.01；2018年相應(yīng)的數(shù)據(jù)與MODIS標(biāo)準(zhǔn)溫度的平均誤差為0.619，方差為2.34，說明本文采用的融合技術(shù)方法獲取的年平均地溫數(shù)據(jù)的準(zhǔn)確度較高。

4 結(jié)論與建議

本文針對青藏高原凍土的特點，提出了一種利用多源遙感數(shù)據(jù)進(jìn)行地表溫度提取的反演方法。該技術(shù)方法通過使用MODIS 8天溫度合成數(shù)據(jù)與Landsat8 反演出的地表溫度進(jìn)行擬合，最終得出地表的年平均溫度；再根據(jù)TTOP模型計算凍融指數(shù)，獲得以Landsat8為基準(zhǔn)的凍土分布，其空間精細(xì)程度與準(zhǔn)確度相較于以往的大尺度制圖有了較大程度的提高，研究成果可為相關(guān)領(lǐng)域的應(yīng)用提供基礎(chǔ)數(shù)據(jù)，技術(shù)方法也可為類似研究提供借鑒。

由于地表溫度反演受到多方面因素的影響，其準(zhǔn)確度仍然存在一定的偏差，因此還需要進(jìn)一步考慮地表溫度反演分析的影響因素及其相關(guān)誤差的來源，在此基礎(chǔ)上提高遙感數(shù)據(jù)的提取精度。