廖 瑩，范繼輝，李怡穎，程根偉

（1. 中國科學院、水利部成都山地災(zāi)害與環(huán)境研究所，四川 成都 610041；2. 中國科學院大學，北京 100049）

凍土是指溫度在0 ℃或0 ℃以下并含冰的各種巖石或土，主要分布在對氣候變化敏感的高緯度和高海拔區(qū)域[1]。2007–2016 年，全球變暖導致全球平均多年凍土溫度增加了0.29 ± 0.12 ℃[2]，青藏高原的多年凍土范圍自1970 年以來縮減了約10 萬km2[3-4]。特別是近地表的正負積溫變化將直接影響多年凍土和季節(jié)凍土的消長。凍融指數(shù)是給定時期內(nèi)的溫度累計值[5]，是近年來較為流行的量化凍土季節(jié)性變化的氣候參數(shù)，主要包括季節(jié)凍結(jié)/融化指數(shù)[6]、近似凍結(jié)/融化指數(shù)[7]和年凍結(jié)/融化指數(shù)[8]等。其中，年凍結(jié)/融化指數(shù)是計算一年內(nèi)所有低于/高于0 ℃的日均溫之和[9]，該方法已廣泛運用在北半球[10-11]、加拿大[12]、挪威[13]、蒙古[14]和我國的青藏高原[15-16]、東北[17-18]和北疆[19-20]等地區(qū)。由于凍融指數(shù)及其變化對地表能量交換以及生態(tài)系統(tǒng)過程有重大影響，因此被作為氣候變化的敏感指示器[21]，在氣候變化評估、凍土制圖以及寒區(qū)工程環(huán)境評估等領(lǐng)域應(yīng)用廣泛[22]。例如，Shi 等[8]利用凍融指數(shù)的連續(xù)空間分布數(shù)據(jù)，結(jié)合氣候預測模型[23]估算了北半球的凍土范圍和面積。Peng 等[12]利用空氣融化指數(shù)和環(huán)境參數(shù)模擬了北半球1971-2000 年的活動層厚度。

西藏高原作為青藏高原的主體，平均海拔4 000 m以上，是我國高寒凍土的主要分布區(qū)，多年凍土面積近70 萬km2，約占整個青藏高原多年凍土面積的一半[24]。在氣候變化的背景下，作為“全球氣候變化敏感區(qū)”的西藏高原表現(xiàn)出整體變暖、變濕的特征[25]。而氣溫的升高對高原凍土影響強烈，加之人類活動引起的負面影響，近30 年活動層厚度呈現(xiàn)整體增大趨勢[26-27]。但由于西藏高原的氣象站點極度缺乏，部分站點建成較晚，整體存在分布稀疏、數(shù)據(jù)缺失等問題，因此氣候變化背景下西藏高原凍融指數(shù)變化的研究還較為欠缺，僅有的幾例研究主要聚焦在雅魯藏布江干流或工程沿線區(qū)域[15-16,28]，這些區(qū)域靠近人類活動區(qū)，受人類活動和工程建設(shè)影響強烈，尚不能代表西藏高原整體的凍融特征。因此本文將在ArcGIS 等軟件的支持下，利用西藏高原共34 個氣象站點的日均氣溫數(shù)據(jù)，對1978-2017 年凍融指數(shù)的時空變化特征進行分析，探討其潛在的影響因子，為下一步估算西藏高原凍土范圍和面積、活動層厚度等提供數(shù)據(jù)，并為西藏高原農(nóng)牧業(yè)生產(chǎn)、草地生態(tài)環(huán)境評估、工程建設(shè)和凍害防治等提供參考。

1 材料與方法

1.1 研究區(qū)概況

西藏高原(26°50′ - 36°53′ N，78°25′ - 99°06′ E)位于北半球中緯度地帶，喜馬拉雅山脈及其以北地區(qū)，青藏高原西南邊陲，面積120.22 萬km2，絕大部分區(qū)域(92%以上)位于海拔4 000 m 以上的高寒生態(tài)區(qū)[29](圖1)。西藏高原具有獨特、多樣的氣候特征，大致表現(xiàn)為西北地區(qū)嚴寒干燥，東南地區(qū)溫暖濕潤。2017 年，全區(qū)年均氣溫在-1.1～12.1 ℃，由東南向西北遞減；干濕分明，年均降水量在117.4～929.4 mm，也呈東南向西北遞減的分布規(guī)律(西藏自治區(qū)氣象局http://xz.cma.gov.cn)；植被從東南向西北依次呈現(xiàn)森林、草甸、草原和荒漠的植被景觀，與氣溫和降水的空間分布高度相關(guān)。

圖1 西藏高原氣象站點分布圖Figure 1 Study areas and distribution of national meteorological stations on the Tibetan Plateau

1.2 數(shù)據(jù)來源及預處理

本研究使用的數(shù)據(jù)包括：①1978-2017 年西藏34 個氣象站點日均溫數(shù)據(jù)(距地面2 m 處氣溫)來源于中國地面氣候資料日值數(shù)據(jù)集(http://data.cma.cn)，由于建站較晚，南木林、貢嘎、墨竹工卡、類烏齊、浪卡子、洛隆和八宿7 站點僅有1992-2017 年的數(shù)據(jù)；②數(shù)字高程模型(digital elevation model, DEM)數(shù)據(jù)(分辨率90 m × 90 m)和植被歸一化指數(shù)(normalized difference vegetation index, NDVI)數(shù)據(jù)(2008-2017)來源于中國科學院資源環(huán)境科學數(shù)據(jù)中心(http://www.resdc.cn)；③人口、地區(qū)生產(chǎn)總值(gross domestic product, GDP)和第一產(chǎn)業(yè)增加值數(shù)據(jù)為34 個氣象站點所在縣(市) 2013 - 2016 年的多年均值，來源于《中國縣域統(tǒng)計年鑒(縣市卷-2014)》、《中國縣域統(tǒng)計年鑒(縣市卷-2015)》、《中國縣域統(tǒng)計年鑒(縣市卷-2016)》和《中國縣域統(tǒng)計年鑒(縣市卷-2017)》[國家數(shù)字電子圖書館(http://www.nlc.cn)]。

考慮到氣象數(shù)據(jù)中存在缺測、漏測等情況，利用缺測站點與周邊站點間的數(shù)據(jù)，按照空間相關(guān)性插值，對缺測部分進行插補。對不同的缺測情況采用不同的插補方法，具體為：①單日缺測：取缺測日前后1 d 平均值；②連續(xù)2 日缺測：第1 缺測日取其前2 d 平均值，第2 缺測日取其后2 d 平均值；③缺測連續(xù)2 日但不超1 年：取前后各1 年同時段的平均值；④缺測日數(shù)超1 年：利用有數(shù)據(jù)年份的數(shù)據(jù)與周邊站點建立相關(guān)關(guān)系，采用最相關(guān)的2～3 個站點數(shù)據(jù)進行插值。

1.3 研究方法

1.3.1 凍融指數(shù)計算

凍融指數(shù)包括凍結(jié)指數(shù)和融化指數(shù)，是凍結(jié)/融化期內(nèi)所有氣溫小于/大于0 ℃的累計值?？紤]到西藏高原冷暖季分明，為保證計算的冷季連續(xù)，本研究以每年7 月1 日至翌年6 月30 日為周期計算凍結(jié)指數(shù)，以一個日歷年的1 月1 日至12 月31 日為周期計算融化指數(shù)[10]。其計算公式分別為：式中：FI(freezing index)、TI(thawing index)分別為凍結(jié)指數(shù)和融化指數(shù)(℃·d)；Ti和Tj分別為負值逐日溫度(℃)和正值逐日溫度(℃)；NF和NT分別為年內(nèi)溫度低于和高于0 ℃天數(shù)。

1.3.2 趨勢分析

采用一元線性回歸分析計算34 個氣象站點氣溫、凍結(jié)指數(shù)和融化指數(shù)的變化趨勢并通過最小二乘法擬合的斜率反映，當斜率 ＞ 0 時，表示增加趨勢，當斜率 ＜ 0 時，表示降低趨勢。斜率計算方法為[30]：

式中：SlopeX為站點x的變化斜率；i為序列序號；n為序列長度；Xi為站點x第i年的值。

1.3.3 空間插值

通過反距離插值法(inverse distance weighted)對34 個氣象站點的凍融指數(shù)進行空間插值[31-33]，以獲得空間連續(xù)的凍融指數(shù)均值和變化趨勢。反距離插值法是對所有樣本點線性加權(quán)，以估計未知點的值。權(quán)重與未知點與樣本點之間的距離成反比，距離越近，權(quán)重越大。具體公式為：

本研究中1978-1991 年的凍融指數(shù)空間結(jié)果由基于拉薩、察隅等27 個站點數(shù)據(jù)空間插值結(jié)果計算得到，1992-2017 年的結(jié)果由基于增加站點后的34 個站點數(shù)據(jù)空間插值計算得到。

1.3.4 相關(guān)關(guān)系分析

采用Pearson 相關(guān)系數(shù)反映凍融指數(shù)與緯度、海拔、人口和GDP 等因子間的相關(guān)程度[34]。計算公式如下：

式中：rij為凍結(jié)(融化)指數(shù)i和某影響因子j之間的相關(guān)系數(shù)，n為樣本數(shù)。

由于各因子間會相互影響，造成相關(guān)關(guān)系結(jié)果的不準確，因此本文利用偏相關(guān)分析(partial correlation analysis)，在控制其他影響因子的基礎(chǔ)上，計算某一影響因子與凍結(jié)(融化)指數(shù)的偏相關(guān)關(guān)系。公式如下[35]：

式中：rij,k為控制某影響因子k的作用后，凍結(jié)(融化)指數(shù)i和另一影響因子j之間的偏相關(guān)系數(shù)?？刂谱兞吭黾訒r的公式以此類推。

2 結(jié)果與分析

2.1 西藏高原氣溫時空變化特征

基于1978-2017 年34 個氣象站點數(shù)據(jù)，計算西藏高原年平均氣溫，并分析其時空變化特征(圖2)。圖2a 表明，以34 個氣象站點代表的西藏高原1978-2017 年平均氣溫在3.04～5.72 ℃，均值4.41 ℃，并呈波動上升趨勢，氣候傾向率為0.06 ℃·a-1。圖2b顯示，近40 年西藏高原年平均氣溫整體呈增加趨勢，類烏齊附近區(qū)域的氣溫增速最快，僅有八宿附近的氣溫呈現(xiàn)下降趨勢。

圖2 1978?2017 年西藏高原年均氣溫年際變化趨勢(a)和空間變化趨勢(b)Figure 2 Inter-annual variation (a) and spatial variation trends (b) of air temperature on the Tibetan Plateau during 1978?2017

2.2 凍融指數(shù)年際間變化特征

西藏高原凍結(jié)指數(shù)的年際變化過程顯示，近40 年西藏高原凍結(jié)指數(shù)介于386.8～886.1 ℃·d，最低值出現(xiàn)在2017 年(386.8 ℃·d)，最高值出現(xiàn)在1982 年(886.1 ℃·d)，均值為591.2 ℃·d (圖3a)；凍結(jié)指數(shù)總體呈現(xiàn)波動下降的趨勢(P＜ 0.01)，下降速率9.2 ℃·d·a-1。凍結(jié)指數(shù)的距平統(tǒng)計結(jié)果(圖3b)表明西藏高原凍結(jié)指數(shù)年際變化十分明顯，凍結(jié)指數(shù)距平值在近40 年內(nèi)逐漸降低(P＜ 0.01)，21 世紀以后，西藏高原各年的凍結(jié)指數(shù)小于近40 年均值。

圖3 1978?2017 年西藏高原凍結(jié)指數(shù)年際變化趨勢(a)和距平變化(b)Figure 3 Inter-annual variation (a) and variations in the departures (b) of the freezing index on the Tibetan Plateau during 1978 ? 2017

由西藏高原融化指數(shù)年際變化過程可以看出，近40 年西藏融化指數(shù)范圍介于1 960.3～2 521.5℃·d，最低值出現(xiàn)在1986 年(1 960.3 ℃·d)，最高值出現(xiàn)在2009 年(2 521.5 ℃·d)，均值為2 209.3 ℃·d(圖4a)。融化指數(shù)整體呈現(xiàn)波動上升的趨勢(P＜0.01)，上升速率為11.1 ℃·d·a-1。融化指數(shù)距平統(tǒng)計圖4b 表明，西藏高原融化指數(shù)整體年際差異較大，融化指數(shù)距平值在近40 年內(nèi)逐漸上升(P＜ 0.01)，在1995-2002 年內(nèi)存在數(shù)值的小幅震蕩，1995 年以前，西藏高原各年的融化指數(shù)小于近40 年平均值，自2002 年以后則大于近40 年平均值。

圖4 1978?2017 年西藏高原融化指數(shù)年際變化趨勢(a)和距平變化(b)Figure 4 Inter-annual variation (a) and variations in the departures (b) of the thawing index on the Tibetan Plateau during 1978 ? 2017

2.3 凍融指數(shù)空間變化特征

西藏高原34 個氣象站點凍結(jié)指數(shù)與融化指數(shù)的統(tǒng)計特征(表1)顯示安多的凍結(jié)指數(shù)多年均值最大(1 759.95 ℃·d)，察隅的凍結(jié)指數(shù)多年均值最小(0.03 ℃·d)。凍結(jié)指數(shù)均值的空間分布特征呈現(xiàn)自東南向西北逐漸遞增的趨勢(圖5a)。1978-2017 年間，除八宿和察隅之外，所有站點的凍結(jié)指數(shù)均呈下降趨勢(表1)，其中，類烏齊的凍結(jié)指數(shù)下降速率最快，傾向率為-17.4 ℃·d·a-1，林芝的凍結(jié)指數(shù)下降速率最慢，傾向率為-0.5 ℃·d·a-1。西藏高原凍結(jié)指數(shù)的傾向率整體呈現(xiàn)自東南向西北遞減的空間分布特征(圖5b)，表明西藏高原西北地區(qū)凍結(jié)指數(shù)下降速率較西藏高原東南地區(qū)更快。

圖5 1978?2017 年西藏高原凍結(jié)指數(shù)均值(a)和變化趨勢(b)Figure 5 Mean annual values (a) and spatial variation trends (b) of the freezing index (FI) on the Tibetan Plateau during 1978 ? 2017

據(jù)西藏高原34 個氣象站點融化指數(shù)統(tǒng)計特征(表1)顯示，察隅的融化指數(shù)多年均值最大(4 446.66℃·d)，安多的融化指數(shù)多年均值最小(940.71 ℃·d)。融化指數(shù)均值在西藏高原東部和中南部個別區(qū)域出現(xiàn)高值中心(圖6a)。1978-2017 年，除八宿附近區(qū)域的融化指數(shù)呈下降趨勢，西藏高原整體經(jīng)歷了融化指數(shù)的快速上升，類烏齊和拉薩的融化指數(shù)上升最快(圖6b)，傾向率分別為22.2 和18.4 ℃·d·a-1(表1)。

圖6 1978?2017 年西藏高原融化指數(shù)均值(a)和變化趨勢(b)Figure 6 Mean annual values (a) and spatial variation trends (b) of the thawing index (TI) on the Tibetan Plateau during 1978 ? 2017

表1 西藏高原34 個氣象站點凍結(jié)指數(shù)與融化指數(shù)的統(tǒng)計特征Table 1 Summary of mean annual freezing and thawing indices for 34 meteorological stations on the Tibetan Plateau

2.4 凍融指數(shù)的影響因子分析

采用Pearson 相關(guān)系數(shù)[36]來衡量凍融指數(shù)與自然因子(如緯度、經(jīng)度、海拔)、NDVI 和社會經(jīng)濟因子(如人口、GDP、第一產(chǎn)業(yè)增加值)共7 個因子間的線性相關(guān)性強弱(表2)。表中顯示，凍結(jié)指數(shù)和融化指數(shù)與自然因子的相關(guān)關(guān)系更強。凍結(jié)指數(shù)與緯度、經(jīng)度、海拔和NDVI 都顯著相關(guān)，與緯度和海拔呈極顯著正相關(guān)關(guān)系(P＜ 0.01)，與經(jīng)度和NDVI 呈顯著負相關(guān)關(guān)系(P＜ 0.05)；融化指數(shù)與經(jīng)度和海拔分別呈現(xiàn)顯著正相關(guān)(P＜ 0.05)和極顯著負相關(guān)關(guān)系(P＜ 0.01)，與緯度和NDVI 關(guān)系不顯著(P＞ 0.05)。凍結(jié)指數(shù)和融化指數(shù)與社會因子之間并無顯著的相關(guān)關(guān)系。Pearson 相關(guān)分析表明，凍結(jié)指數(shù)與海拔呈現(xiàn)最高正相關(guān)關(guān)系，融化指數(shù)與海拔呈現(xiàn)最高負相關(guān)關(guān)系。

凍融指數(shù)與7 個因子的偏相關(guān)結(jié)果(表2)表明凍結(jié)指數(shù)與緯度和海拔呈極顯著正相關(guān)關(guān)系(P＜0.01)；融化指數(shù)與海拔呈極顯著負相關(guān)關(guān)系(P＜0.01)，與NDVI 呈顯著負相關(guān)(P＜ 0.05)。凍結(jié)指數(shù)與海拔呈現(xiàn)最高正相關(guān)關(guān)系，融化指數(shù)與海拔呈現(xiàn)最高負相關(guān)關(guān)系。

表2 西藏高原凍融指數(shù)與影響因子的Pearson 相關(guān)及偏相關(guān)系數(shù)Table 2 Pearson correlation and partial correlation coefficients of the freezing and thawing indices and influencing factors on the Tibetan Plateau

綜合以上兩種分析結(jié)果，認為凍融指數(shù)與海拔的相關(guān)性最強。將凍結(jié)指數(shù)和融化指數(shù)與海拔的偏相關(guān)關(guān)系以散點圖的形式更直觀地表達(圖7)，可以觀察到凍結(jié)指數(shù)與海拔有很強的正線性關(guān)系(P＜0.01) (圖7a)，融化指數(shù)與海拔具有很強的負線性關(guān)系(P＜ 0.01) (圖7b)，R2值均超過了0.5。

圖7 西藏高原凍結(jié)指數(shù)(a)和融化指數(shù)(b)與海拔的偏相關(guān)散點圖Figure 7 Scatterplot of the partial correlation between the freezing index (a) and thawing index (b) and altitude on the Tibetan Plateau

3 討論

在全球變化背景下，青藏高原地區(qū)氣溫持續(xù)上升已是得到觀測驗證的事實[37-38]，本研究中計算的西藏高原氣溫傾向率0.06 ℃·a-1與SROCC 評估[39]表明的增溫速率0.03 ℃·a-1和楊春艷等[40]得到的西藏近50 年氣溫傾向率每10 年升高0.58 ℃的結(jié)果基本一致。西藏高原1978-2017 年凍融指數(shù)及其變化趨勢均符合環(huán)北極地區(qū)和基于氣象站點觀測數(shù)據(jù)的青藏高原凍融指數(shù)范圍和傾向率(表3)，由于不同數(shù)據(jù)來源、區(qū)域范圍差異和所選氣象站點不同，存在具體數(shù)值上的偏差。值得注意的是，基于不同數(shù)據(jù)來源的青藏高原凍融指數(shù)及其變化趨勢差異較大，可能是因為全球格網(wǎng)氣候數(shù)據(jù)產(chǎn)品的空間分辨率較大(2.5° × 2.5°/0.5° × 0.5°)[8]，更易忽視下墊面的復雜程度，尤其是青藏高原地區(qū)地勢起伏劇烈，低分辨率的數(shù)據(jù)產(chǎn)品難以模擬復雜的氣溫分布情況，導致計算結(jié)果的出入較大。因此，針對尺度較小的區(qū)域，利用氣象站點的觀測數(shù)據(jù)計算凍融指數(shù)會更接近真實值。

表3 不同數(shù)據(jù)來源和不同區(qū)域的凍融指數(shù)及其變化趨勢對比Table 3 Comparison of freezing and thawing indices and the variation trends of different data sources and different regions

本研究基于各站點的凍融指數(shù)及其變化趨勢進行空間插值以獲得區(qū)域內(nèi)凍融指數(shù)連續(xù)空間分布特征。凍結(jié)指數(shù)均值呈現(xiàn)由東南向西北逐漸增加的特征與氣溫的分布特征相關(guān)，而對氣候變化的敏感度更高也是導致藏北地區(qū)凍結(jié)指數(shù)下降速率較西藏高原其他地區(qū)更高[40]的主要原因。融化指數(shù)在西藏絕大部分區(qū)域都呈上升趨勢，并在拉薩、林芝和類烏齊等區(qū)域呈現(xiàn)出點狀高值中心，可能與近年來農(nóng)業(yè)開墾、城市擴張、產(chǎn)業(yè)發(fā)展等較強的社會經(jīng)濟活動有關(guān)，同時，這些區(qū)域也是西藏生態(tài)退化區(qū)[29]，對人類活動的影響更為敏感。

本研究歸納的4 個自然因子與西藏高原凍融指數(shù)都有一定的相關(guān)關(guān)系，但納入考慮的3 個社會經(jīng)濟因子與凍融指數(shù)并無較顯著的相關(guān)關(guān)系，說明凍融指數(shù)作為一項氣候參數(shù)更易受到自然條件的影響，尤其是與氣溫關(guān)系緊密的自然因子。由于緯度位置決定地球上的熱量分配，近地面氣溫也存在垂直高度上的變化規(guī)律，所以緯度和海拔與氣溫的關(guān)系最緊密，并與凍融指數(shù)顯著相關(guān)。由于西藏高原氣象站點極為稀少，僅有的34 個標準氣象站也大多分布在海拔3 500～4 500 m 的位置，因此本研究利用空間插值獲得的凍融指數(shù)連續(xù)空間分布特征所能代表的范圍和精度有限。雖考慮到凍融指數(shù)與海拔具有極顯著相關(guān)關(guān)系，但目前成熟的插值方法尚未能將各影響因素的貢獻度作為插值過程中的權(quán)重，下一步研究將嘗試插值方法上的創(chuàng)新，優(yōu)化凍融指數(shù)及其變化趨勢的空間結(jié)果。

土壤的凍脹、融沉及其年內(nèi)循環(huán)會造成土壤力學性質(zhì)的變化，嚴重時會發(fā)生工程凍害，影響土壤穩(wěn)定性和上部建筑的安全[42]。凍融指數(shù)是繪制多年凍土和季節(jié)性凍土分布圖[23]、估算土壤凍結(jié)深度[43-44]、測量凍土活動層厚度[45-46]以及指示凍融作用強度和特征的重要指標，近年來在寒區(qū)工程環(huán)境評估領(lǐng)域應(yīng)用廣泛[47]。隨著西藏自治區(qū)人口增長、工農(nóng)業(yè)發(fā)展、旅游開發(fā)、公路鐵路網(wǎng)逐步完善與基礎(chǔ)設(shè)施建設(shè)的開展，預防工程凍害是重中之重，本研究依據(jù)最新的氣象資料計算分析凍融指數(shù)并獲得連續(xù)的空間分布數(shù)據(jù)，可以為下一步估算西藏高原多年凍土面積、模擬多年凍土空間分布范圍、測量各區(qū)域土壤凍結(jié)深度等提供重要參數(shù)數(shù)據(jù)，為西藏高原農(nóng)牧業(yè)生產(chǎn)、草地生態(tài)環(huán)境評估、工程布局和凍害治理等提供科學參考。

4 結(jié)論

通過上述分析，可以得出以下幾點主要結(jié)論：

1) 1978-2017 年西藏高原的氣溫在波動中上升，多年平均值為4.41 ℃，氣候傾向率為0.06 ℃·a-1，約是全球同期升溫率的2 倍，對全球變化高度敏感。

2) 1978-2017 年西藏高原凍結(jié)指數(shù)范圍在386.8～886.1 ℃·d，平均值為591.2 ℃·d，融化指數(shù)范圍在1 960.3～2 521.5 ℃·d，平均值為2 209.3 ℃·d。凍結(jié)指數(shù)和融化指數(shù)的變化趨勢相反，分別隨時間波動下降和波動上升，變化速率為-9.2 和11.1 ℃·d·a-1。

3) 近40 年凍結(jié)指數(shù)均值和下降速率均呈現(xiàn)出自東南向西北遞增的空間分布規(guī)律。融化指數(shù)均值較高和上升明顯的區(qū)域集中在以拉薩為主的西藏高原中南部和東部。

4) 海拔是影響西藏高原凍融指數(shù)的關(guān)鍵因子，與凍結(jié)指數(shù)呈現(xiàn)極顯著正相關(guān)關(guān)系，與融化指數(shù)呈現(xiàn)極顯著負相關(guān)關(guān)系。