1978?2017 年西藏高原凍融指數時空變化特征

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

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

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

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

1 材料與方法

1.1 研究區概況

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

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

1.2 數據來源及預處理

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

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

1.3 研究方法

1.3.1 凍融指數計算

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

1.3.2 趨勢分析

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

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

1.3.3 空間插值

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

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

1.3.4 相關關系分析

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

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

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

式中：rij,k為控制某影響因子k的作用后，凍結(融化)指數i和另一影響因子j之間的偏相關系數。控制變量增加時的公式以此類推。

2 結果與分析

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

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

圖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 凍融指數年際間變化特征

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

圖3 1978?2017 年西藏高原凍結指數年際變化趨勢(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

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

圖4 1978?2017 年西藏高原融化指數年際變化趨勢(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 凍融指數空間變化特征

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

圖5 1978?2017 年西藏高原凍結指數均值(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

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

圖6 1978?2017 年西藏高原融化指數均值(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 個氣象站點凍結指數與融化指數的統計特征Table 1 Summary of mean annual freezing and thawing indices for 34 meteorological stations on the Tibetan Plateau

2.4 凍融指數的影響因子分析

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

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

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

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

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

3 討論

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

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

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

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

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

4 結論

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

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

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

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

4) 海拔是影響西藏高原凍融指數的關鍵因子，與凍結指數呈現極顯著正相關關系，與融化指數呈現極顯著負相關關系。