基于多源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為基準的凍土分布，其空間精細程度與準確度相較于以往的大尺度制圖有了較大程度的提高，研究成果可為相關領域的應用提供基礎數據，技術方法也可為類似研究提供借鑒。

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