1901—2018年三江源地區大氣凍融指數變化特征

陳方方， 羅棟梁， 劉 磊， 金會軍，3， 李 超

（1.中國科學院西北生態環境資源研究院凍土工程國家重點實驗室，甘肅 蘭州 730000； 2.中國科學院大學，北京 100049；3.東北林業大學土木工程學院寒區科學與工程研究院，黑龍江 哈爾濱 150040； 4.蘭州交通大學測繪與地理信息學院，甘肅 蘭州 730070）

0 引言

凍土是高寒生態系統的重要組成部分，其變化深刻影響地表水熱平衡、生物地球化學循環乃至整個高寒生態系統的穩定［1］。在氣候變暖和人類活動增強疊加作用下，凍土發生了顯著退化，表現為地溫升高、季節凍深變淺、活動層增厚、融區貫通或擴展、多年凍土厚度減薄乃至消失［2-5］，由此引起生態系統穩定性減弱和資源環境承載壓力增大等問題［6］。凍結融化指數是凍土研究相關的重要參數之一，也是判斷凍土存在狀態與長期變化的重要指標之一［7］。在季節凍土區，年凍結指數是計算年最大凍結深度的重要依據，以評估寒區工程施工的凍脹量等［8］。在多年凍土區特別是資料匱乏的地區，年凍融指數是計算多年凍土和活動層厚度時空分布變化的關鍵因子之一［9］。

凍土區因交通不便和氣候嚴酷造成實測地氣溫資料匱乏，導致凍融指數研究以前多基于國家基準臺站或自動氣象站的觀測資料計算。趙紅巖等［10］基于地面溫度觀測資料計算了青藏鐵路沿線凍融指數，其分析表明地面凍融指數主要受緯度和海拔的影響。姜逢清等［11］研究了青藏鐵路沿線凍融指數變化趨勢，結果表明1966—2004年凍結指數減小，而融化指數增加。Wu 等［12］研究發現青藏高原多年凍土區的凍結指數減少速率高于季節凍土區，而融化指數增加速率則小于季節凍土區。劉磊和羅棟梁［13］基于國家基準臺站數據分析了雅魯藏布江流域大氣/地面凍融指數時空變化，結果表明1977—2017年大氣/地面凍結指數均減小，融化指數均上升。Luo等［14］計算了我國東北北部1972—2005年年平均氣溫和地面溫度、大氣/地面凍融指數，結果顯示在過去40年氣候顯著變暖背景下，凍結指數減小，而融化指數顯著增加。綜上可知，以往研究大多基于國家基準臺站觀測資料，其空間分布極不均勻，計算時段也偏短，針對整個研究區長時間序列的凍融指數變化研究較少；且大都集中于資料較豐富地區，難以充分反映某一地區凍融指數的整體時空變化情況。

三江源地區位于青藏高原中東部，是長江、黃河、瀾滄江的發源地，其特殊的地理位置、豐富的自然資源、重要的生態服務功能使其成為我國重要生態安全屏障。目前三江源地區的植被和水源涵養等高寒生態水文環境要素及經濟社會可持續發展不僅得到地方政府和行業部門重視，也引起相關科學研究工作者的廣泛關注［15-17］。不少研究指出，高寒植被、生物生產力和土壤化學元素循環等相關生態環境要素的變化與凍土環境及其變化有著緊密聯系［18］。因此，研究三江源地區凍融指數的時空變化特征，對于了解該區凍土變化及其與生態環境的協同變化具有較重要意義。本文基于英國東英吉利大學（University of East Anglia）大學氣候研究中心（Climatic Research Unit，CRU）的逐月氣溫數據集，并利用國家基準臺站數據驗證，計算并分析了三江源地區大氣凍融指數的時空變化特征，研究結果可在一定程度上反映其氣候和凍土變化特征，同時可為該地區高寒生態環境長時期變化研究提供科學依據。

1 數據和方法

1.1 研究區概況

三江源地區包括長江源區、瀾滄江源區和黃河源區，是我國重要的固碳、水源涵養地和生態屏障區［19］（圖1），其中長江源區以直門達為流域出口、瀾滄江源區以昌都為流域出口、黃河源區以唐乃亥為流域出口［20］。區內凍土、冰川、湖泊、高寒濕地等廣泛分布，是長江、黃河和瀾滄江流域重要的水源地和補給區，被譽為“中華水塔”［21］。區域氣候為典型高原大陸性氣候，冬季嚴寒漫長，夏季溫暖短暫；溫度日較差大而年較差小，雨熱同季，干濕季分明［22］。年均氣溫-5.6～7.8 ℃，年降水量262.2～772.8 mm，年日照時數2 300～2 900 h［6］。因位于大片連續多年凍土區邊緣，三江源地區大片連續、不連續及島狀多年凍土和季節凍土交替分布。

圖1 三江源地區地理范圍及國家基準臺站分布Fig.1 Geographical background and the national meteorological stations distributed in the Three-River Source Region

1.2 資料來源

由于該區氣象站稀少且分布不均，且大多在1950 年后才建站觀測，因而不便利用臺站觀測數據開展凍融指數長時間序列的時空變化研究。東英吉利大學（University of East Anglia）氣候研究中心（Climatic Research Unit）CRU TS 4.03（以下簡稱CRU）通過整合已有的若干知名數據集，重建了一套覆蓋完整、高分辨率且無缺測的月平均地表氣候要素數據集，成為一套較為完整的全球百年尺度氣候變化數據集。這套資料和中國已有氣候數據相比具有如下優點：第一，青藏高原20 世紀前半期器測資料嚴重缺乏，CRU 資料盡管因復雜地形和插值方法本身問題存在誤差，經比較仍可作為有一定信度的參考資料；第二，中國現有百年溫度序列的時間分辨率為年或季，而CRU 資料時間分辨率為月，時間分辨率較高；第三，這個溫度數據集基于觀測結果統計內插得到，故減少了代用資料的不確定性［23］。本文所選取的CRU 逐月氣溫空間分辨率為0.5°×0.5°，時間范圍為1901年1月至2018年12月；共使用174個格點的逐月氣溫資料分析三江源地區凍結指數和融化指數的時空變化。

為檢驗CRU 資料在三江源地區的適用性，本研究從中國氣象數據網（http：//data.cma.cn）下載了中國地面氣候資料月值數據集進行驗證。三江源地區共計19 個國家基準臺站，考慮到部分臺站被撤銷、遷移及觀測時間不連續等問題，總共選取了12個站點1960 年1 月至2018 年12 月的氣溫數據分析驗證（表1）。評價方法是將國家基準臺站的逐月氣溫與對應CRU 格點1960—2018 年的逐月氣溫做統計分析。由于CRU空間分辨率為0.5°，單個格點代表空間范圍較大，因此考慮海拔差異和氣溫垂直遞減率賦予CRU 數據月均溫相應增量［-0.402 ℃·（100m）-1］［24］，在此基礎上比較基準臺站和CRU 格點逐月氣溫統計差異。具體賦值過程為以基準臺站海拔為基準，根據其與對應CRU 格點的海拔差值，賦予CRU逐月氣溫相應增量。

表1 三江源地區國家基準臺站與CRU格點逐月氣溫差值Table 1 Monthly air temperature difference between CRU grid and the national meteorological stations in the Three-River Source Region

1.3 研究方法

1.3.1 凍融指數計算

大氣凍結/融化指數（以下均簡稱凍結/融化指數）的定義為：某個給定時段內日平均氣溫低于/高于0 ℃的持續時間和其數值乘積的總和，以℃·d或℃·mon 表示。本研究基于逐月氣溫計算凍融指數，其計算公式如下：

式中：FI、TI分別為凍結指數和融化指數（℃·d）；Ti為第i月的月平均氣溫（℃）；Di為第i月天數。本研究中，根據每個格點1901年1月到2018年12月的逐月氣溫計算每個格點1901—2018 年的年凍融指數。根據每個格點所有年份的凍融指數求其平均值，得到1901—2018年的三江源地區平均凍融指數。

1.3.2 滑動平均法

滑動平均法又稱移動平均法，是對時間序列數據進行平滑和濾波的一種處理方法。通過計算一定滑動長度的平均值，目的是消除統計序列中的隨機波動，用平均值顯示時間序列的變化趨勢［25］。本文利用滑動平均法對凍結/融化指數1901—2018 年時間序列數據進行平滑和濾波。對于樣本量為n的序列y，其滑動平均序列可表示為：

式中：k為滑動長度。本文利用Python 3.7編程實現三江源地區1901—2018 年凍融指數的滑動平均處理。根據1901—2018 年所有格點凍結/融化指數平均值，得到三江源地區1901—2018 年凍結/融化指數年際變化。同時選擇k=11 年對三江源地區凍結/融化指數的年際變化做滑動平均分析，并根據分析結果對凍結/融化指數滑動平均曲線分段線性擬合。

1.3.3 線性傾向估計

利用一次線性傾向分析凍結/融化指數的變化特征。其公式為：

式中：y為凍結/融化指數；x為年份；a為回歸常數；b為反映凍結/融化指數變化速率的回歸系數（顯著水平P＜0.05）。本文利用Python語言的numpy程序包對三江源地區1901—2018 年凍結/融化指數進行一次線性傾向估計，獲得每個格點凍結/融化指數變化的線性傾向率，進而通過空間分析獲得整個三江源地區凍融指數線性變化的區域差異。考慮到凍結/融化指數自身的空間差異性，本文通過相對變化量來反映每個格點凍融指數相對其自身的變幅，將1901—2018 年的差值與平均值的比值作為相對變化量。

1.3.4 凍結融化指數空間分布

基于各格點凍結/融化指數平均值與其地理信息的多元統計關系，利用ArcMap 空間插值獲得其空間分布特征。對每個格點1901—2018 年的凍結/融化指數求平均值，得到其多年平均值。利用三江源地區500 m 分辨率高程數據重采樣獲得對應格點平均海拔。基于每個格點的經緯度、海拔與凍結融化指數建立其與經緯度、海拔的多元統計關系。同時，為消除量綱影響，對海拔數值乘以0.01 之后再進行回歸分析。最終獲得如下關系：

式中：FI、TI分別為凍結指數和融化指數；x、y、h分別為經度、緯度和海拔。在ArcMap 中基于多元線性關系對凍結和融化指數進行空間插值。

2 結果與分析

2.1 CRU數據與國家基準臺站數據的比較

CRU 數據集已在北半球及中國部分地區的氣候變化和多年凍土研究中廣泛使用，具有較高的精度和適用性［26-27］。在青藏高原地區，CRU 資料也有一定的可信度。本文利用相關系數（CRU 數據和國家基準臺站數據相關程度）、解釋方差R2（CRU 數據可解釋國家基準臺站數據的方差百分比）、均方根誤差RMSE（CRU 數據與國家基準臺站數據偏差的平方和與觀測次數比值的平方根）及平均絕對誤差MAE 來分析CRU 數據在三江源地區的適用性。結果表明，國家基準臺站的逐月氣溫與對應CRU 格點逐月氣溫數據具有較好的一致性（圖2）。兩者變化范圍均為-30～20 ℃，RMSE=3.29，MAE=1.89。相關系數0.92 及解釋方差R2=0.93 表明CRU 逐月氣溫與國家基準臺站數據具有很強的正相關性，適用于三江源地區氣溫及凍結融化指數的時空變化分析。

圖2 CRU逐月氣溫與國家基準臺站逐月氣溫相關分析Fig.2 Correlation of monthly air temperatures between the CRU and the national meteorological stations in the Three-River Source Region

2.2 凍結指數時空變化特征

三江源地區多年平均凍結指數的空間分布較好地體現了其與海拔的相關性：海拔自西向東逐漸遞低，凍結指數也隨之減少（圖3）。同時多年凍土的連續性自西向東降低，由大片連續多年凍土過渡到不連續及島狀多年凍土甚至季節凍土，凍結指數在空間分布上也呈現了明顯的自西向東遞減趨勢：由最西端的3 400 ℃·d減小到最東端的不到600 ℃·d，平均值為1 940 ℃·d，大部分地區的凍結指數值界于800～2 600 ℃·d。其中凍結指數最大值出現于長江源區西部，最小值則分布在黃河源區東南部。

圖3 三江源地區1901—2018年多年平均凍結指數空間分布Fig.3 Spatial distribution of average air freezing index from 1901 to 2018 in the Three-River Source Region

三江源地區1901—2018 年凍結指數整體呈下降趨勢，其遞減斜率為-1.1 ℃·d·a-1，但并非呈一致性下降趨勢，而是隨氣候波動經歷了下降-上升-下降等3 個波動變化階段（圖4）。其中，1901—1943年為波動下降階段，其變化趨勢為-3.4 ℃·d·a-1；1943—1966 年呈波動上升趨勢，其變化趨勢為8.8 ℃·d·a-1；1966—2018 年又呈波動下降趨勢，其變化趨勢為-4.3 ℃·d·a-1。

圖4 三江源地區1901—2018年凍結指數時間序列變化Fig.4 Time series of air freezing index from 1901 to 2018 in the Three-River Source Region

三江源地區1901—2018 年凍結指數線性傾向率變化范圍-2.6～-0.11 ℃·d·a-1［圖5（a）］，在空間上自西向東呈減少趨勢，以長江源區的變化速率最大，黃河源區的變化速率最小。除長江和黃河源區北部凍結指數變化不顯著外，其他區域均顯著變化。由圖5（b）可知，凍結指數相對變化量具有自南向北減少趨勢，黃河源區凍結指數的相對變化量較小。總體而言，三江源地區凍結指數的相對變化量均較小，最小為-0.1，表明三江源地區凍結指數的變化幅度不大。

圖5 三江源地區1901—2018年凍結指數線性傾向（a）（實心圓表示顯著，P＜0.05）及其相對變化量（b）區域差異Fig.5 Spatial difference of the changing rates（a）（solid circle denotes the change is significant，P＜0.05）and the relative change（b）of the air freezing index from 1901 to 2018 in the Three-River Source Region

2.3 融化指數時空變化特征

融化指數的空間分布與凍結指數呈相反趨勢，即隨海拔升高而減小，隨海拔降低而增大，故三江源地區西部融化指數較小，而東部融化指數較大，在空間分布上自西向東遞增（圖6）。三江源地區融化指數變化范圍為0～1 800 ℃·d，平均值為787 ℃·d。其中長江源區融化指數最小，而黃河源區融化指數最大。

圖6 三江源地區1901—2018年多年平均融化指數空間分布Fig.6 Spatial distribution of average air thawing index from 1901 to 2018 in the Three-River Source Region

三江源地區融化指數在1901—2018 年整體呈波動上升趨勢，其遞增斜率為0.34 ℃·d·a-1，具體而言經歷了上升-下降-上升等3 個階段（圖7）。其中，1901—1943 年為波動上升階段（3.3 ℃·d·a-1），而1943—1981 年為波動下降階段（-3.1 ℃·d·a-1），1981—2018 年為波動上升階段（2.9 ℃·d·a-1）。融化指數上升和下降的斜率大致相等，表明變化速率大致相同。

圖7 三江源地區1901—2018年融化指數時間序列變化Fig.7 Time series of air thawing indices from 1901 to 2018 in the Three-River Source Region

三江源融化指數1901—2018 年間的線性傾向率變化范圍為-0.29～0.95 ℃·d·a-1，且大部分地區融化指數變化斜率大于零，表明融化指數呈上升趨勢，只有小部分地區呈下降趨勢［圖8（a）］。在空間分布上，融化指數的變化率自西向東呈減小趨勢，以長江源區融化指數遞增速率最大，黃河源區的遞增速率最小。融化指數相對變化量的變化范圍較大，為0.07～0.86，大多集中在0.1～0.3［圖8（b）］。從空間分布來看，三江源地區融化指數1901—2018年的相對變化量以西部長江源頭和中部最大，黃河源區東部融化指數的相對變化量較小。

圖8 三江源地區1901—2018年融化指數線性傾向（a）（實心圓圈表示顯著水平P＜0.05）和相對變化量（b）區域差異圖Fig.8 Spatial difference of the change rate（a）（solid circle denotes the change is significant，P＜0.05）and the relative change（b）of the air thawing index from 1901 to 2018 in the Three-River Source Region

3 討論

3.1 CRU數據的適用性

氣候變化研究通常需使用百年以上長時間序列的氣候資料，而青藏高原器測資料時間短、觀測站稀缺且分布極不均勻，因此長時間序列氣候和環境變化以及多年凍土的變化常基于ERA、CRU、MERRA2 等再分析資料進行研究［28-30］。目前CRU逐月氣溫再分析資料在凍土變化研究中仍存在不足。一方面，CRU 逐月氣溫是月尺度資料，相較于日均氣溫在凍融指數計算方面存在一定偏差，特別是凍結和融化過程頻繁交替的夏初或秋冬月份偏差尤其大。當月內逐日氣溫均高于或低于0 ℃時，基于逐月氣溫的計算自然可真實反映凍融指數變化；但在凍融過程頻繁交替時，即逐日氣溫在月內圍繞0 ℃反復波動時，基于逐月氣溫則較難獲得真實凍融指數。青藏高原作為典型的高海拔凍土區，日內溫差變化大，其在春季和秋季的凍融交替過程中，溫度在0 ℃附近頻繁波動，因此基于CRU 等再分析資料獲取準確的凍融指數比高緯度地區面臨更嚴峻的挑戰［31］。另一方面，由于CRU 數據空間分辨率為0.5°×0.5°，每個格點所包含的實地面積約為3 000 km2，而青藏高原地形異常復雜，在較小空間范圍內海拔變化可達數百米乃至上千米。考慮到氣溫的垂直遞減率，基于再分析資料計算凍融指數勢必存在低估或高估的現象；與青藏高原相比，CRU 數據在北半球其他地區的適用性更高［31］。若要精準獲得青藏高原地區凍融指數的時空變化規律，需在考慮復雜地形和局地因子影響的基礎上對現有再分析資料降尺度以發展更高時空分辨率的氣候數據。在今后研究工作中，擬將CRU 逐月氣溫與已有逐日氣溫再分析資料或空間分辨率更高的再分析資料［25，32］與地形因子校正、實地觀測資料相結合，以盡量減小誤差。

3.2 凍融指數的時空變化及其影響

對比凍融指數的空間分布，可以發現凍結和融化指數的空間分布大致呈相反趨勢，即自西向東，隨著海拔逐漸下降及多年凍土連續性逐漸降低，凍結指數遞減，融化指數遞增；表明從長江源區到黃河源區年平均氣溫大致在上升。由凍結指數與經緯度和海拔的多元線性關系可知，凍結指數與緯度和海拔正相關，而與經度負相關（圖1和圖3）。三江源地區海拔自西向東在不斷降低，凍結指數與經度的負相關從一定程度上可理解為凍結指數與海拔正相關。因此，三江源地區凍結指數的空間分布呈現出海拔越高、凍結指數就越大。由圖3 可知，由南向北，凍結指數由小于800 ℃·d逐漸增加到2 600 ℃·d，體現了凍結指數在空間分布上的緯度效應。這與曹斌等［33］在研究黑河流域凍結指數的空間分布，何彬彬等［34］研究北疆地區凍結指數的影響因素的結果一致，即凍結指數主要受緯度和海拔的雙重影響。

凍結指數為一年中的負積溫，融化指數則為一年中的正積溫。凍融指數的變化不僅反映了土壤凍結狀態的變化，也在一定程度上反映了氣溫變化和多年凍土對氣候變化的響應幅度。氣溫變化是導致凍融指數變化的直接原因，其劇烈程度與凍融指數的變化幅度相一致。本文研究表明三江源地區1901—2018 年間凍結指數呈波動下降、融化指數呈波動上升趨勢，這大體與易湘生等［27］、陳曉光等［35］、張士鋒等［20］的結果一致，即三江源地區自二十世紀五六十年代以來氣溫呈波動上升趨勢，以冬季增溫最顯著。本文的研究表明，凍融指數氣候傾向率以長江源區和瀾滄江源區的變化速率最大，大致呈自西向東遞減的趨勢，這表明長江源區和瀾滄江源區是整個三江源地區氣溫升幅最大的地區。這與竇睿音［36］對三江源地區近半個世紀的氣候變化研究相一致，即瀾滄江源區與長江源區增溫速率最快，分別為0.35 ℃·d·a-1、0.32 ℃·d·a-1；而黃河源區增溫速率最小，僅為0.19 ℃·d·a-1。此外，凍融指數的氣候傾向率隨海拔升高而增大，反映出地溫較低、連續性較高的多年凍土其退化更嚴重。Ran等［37］發現青藏高原海拔3 600 m 處年均氣溫增溫速率為0.33 ℃·d·a-1，而海拔5 200 m 處的增溫速率為0.49 ℃·d·a-1，這與三江源地區自西向東海拔不斷降低、高海拔大片連續多年凍土區增溫速率最大的結果較為一致。

凍融指數高低反映了土壤的凍融起止時間及持續時長，因此對生態系統、地表水文過程、農業生產、基礎設施等都會產生重要影響［8］。凍融指數作為生態系統的一個重要指示因子，受季節凍融循環影響明顯的生態系統過程可分為生長季、土壤凍結初期、凍結期和融凍期等四個關鍵時期［38］。凍融指數通過決定植物生長發育時期的時長，從而進一步影響動植物的生長進程、動物群落的結構和功能以及相關土壤生化特性。融化指數增大、凍結指數減小，意味著融化期延長、凍結期縮短，表明植被生長季將延長，植物生產力可能得到提高，高寒植被的生態服務功能也將提升；反之，則植被生長季縮短，植物生產力減弱，其生態服務功能亦可能削弱。整體上，三江源地區1901—2018 年間凍結指數減小、融化指數增大，表明高寒植被在過去100 年間可能有了更長的生長季，因而高寒生態環境可能有了一定程度的改善。

4 結論

本文基于國際上廣泛使用的CRU 再分析資料，對1901—2018 年間三江源地區的凍結和融化指數的時空變化及其差異進行了比較分析，主要結論如下：

（1）三江源地區凍結指數的變化范圍為600～3 400 ℃·d，大多數凍結指數集中分布于800～2 600 ℃·d；在空間分布上，緯度和海拔越高，凍結指數就越大。融化指數的變化范圍為0～1 800 ℃·d，空間分布呈現自西向東遞增的趨勢。

（2）三江源地區凍結指數在1901—2018 年整體呈波動減少趨勢，遞減斜率為-1.1 ℃·d·a-1，大致經歷了下降（1901—1943 年）-上升（1943—1966 年）-下降（1966—2018 年）等三個階段；融化指數整體波動增加，遞增斜率為0.34 ℃·d·a-1。

（3）三江源地區凍結指數線性傾向斜率變化范圍-2.6～-0.11 ℃·d·a-1。在1901—2018 年間呈降低趨勢。融化指數線性傾向斜率變化范圍-0.29～0.95 ℃·d·a-1，大部分地區融化指數呈上升趨勢。自西向東，凍結指數遞減斜率在減小，融化指數遞增速率也在減小。凍結指數相對變化量都較小，最小為-0.1，說明三江源地區凍結指數的變化幅度不大。而融化指數相對變化量的變化范圍較大，為0.07～0.86 ℃·d·a-1，但大多集中在0.10～0.30 ℃·d·a-1。

（4）凍融指數的時空變化在一定程度上反映了氣候變化對高寒植物生長的影響，三江源地區凍融指數過去100年的時空變化特征可能意味著高寒植被有了更長的生長季。