中國東北興安嶺地區年凍融頻次的分布規律

魏 豪， 張 澤， Andrey MELNIKOV， 金豆豆， 高思如， 馮文杰

（1.中國科學院西北生態環境資源研究院凍土工程國家重點實驗室，甘肅蘭州730000； 2.東北林業大學土木工程學院/寒區科學與工程研究院，黑龍江哈爾濱150040； 3.東北林業大學東北多年凍土區地質環境系統教育部野外科學觀測研究站/東北多年凍土區環境、道路建設與養護協同創新中心，黑龍江哈爾濱150040； 4.俄羅斯科學院西伯利亞分院梅爾尼科夫凍土研究所，俄羅斯聯邦雅庫茨克117997； 5.中國科學院大學，北京100049）

0 引言

興安嶺地區位于中國東北北部，包括從內蒙古赤峰市到黑龍江漠河市的大興安嶺，以及北起黑龍江岸南至松花江岸的小興安嶺，有中國面積最大、保存相對完好的原始森林，是重要的林業基地和東北邊疆的重要門戶。興安嶺地區包含大部分多年凍土區和部分季節凍土區，與俄羅斯貝加爾湖和遠東地區的不連續多年凍土以及蒙古國山區多年凍土相鄰，是中國高緯度多年凍土的主要分布區和寒區工程建設的重要地區。

近年來全球氣候變暖日趨嚴重，IPCC 報告指出在1880—2012年間，全球的近地層年平均氣溫升高了0.85 ℃［1-2］，中國東北地區的多年凍土和季節凍土都有不同程度的退化，多年凍土區和季節凍土區的南緣都有明顯的北移現象［3］。20 世紀50 年代至21世紀初，東北地區多年凍土總面積從4.8×105km2減少為3.1×105km2，凍土南緣向北移動0.1°～1.1°，平均海拔升高160.5 m［4］。同樣地，東北地區地表凍結指數的下降和融化指數的上升與全球變暖、凍土退化有著相同的趨勢［5-6］。因此，對興安嶺地區的凍土變化及寒區工程建設進行調查研究尤為重要。

俄羅斯凍土學家Ершов［7］在《寒區巖土成因論》中提到，“在自然條件下，一年當中寒區表層的巖土會發生100 余次凍結-融化過程（溫度正負變化100余次）。”在一年內發生凍融循環的次數即為年凍融頻次，在寒區工程建設和使用過程中，凍融循環是對各類材料的力學性質影響最大的因素，因此年凍融頻次的變化對寒區工程的建設、發展以及長期穩定性有著極其重要的影響。

在寒區研究方面，多集中在凍融指數、凍土分布等方面。凍融指數是給定時期內的溫度累計值，是對凍土季節性變化進行量化的氣候參數。Zhang等［8］結合多年凍土區的年平均地表溫度和凍結指數進行研究，發現東北多年凍土區在20 世紀80 年代至90 年代凍土退化嚴重，普遍由低海拔向高海拔、由連續凍土向不連續凍土退化，在此之后多年凍土區面積略有增加；任景全等［9］通過對吉林省季節凍土區凍融指數的研究發現，在空間上吉林省的凍結指數呈現由北向南逐漸減小的趨勢，在時間上凍結指數呈現逐漸減小的趨勢。

在室內試驗方面，多著眼于工程材料，如巖石、混凝土等在多次凍融循環下的強度等力學性質的變化情況。Jamshidi 等［10］、方麗莉等［11］研究發現凍融循環會通過改變土的結構而導致土的工程性質的改變，進而影響基礎設施的耐久性，甚至會導致工程失穩等問題；有學者指出，長期的凍融循環對巖石產生凍融風化、開裂等作用，對許多巖石工程的穩定性產生顯著影響［12-15］。這使得表述凍土分布和寒區氣候變化的凍融指數在工程實踐中存在一定的局限性。而年凍融頻次是定量表述寒區凍融循環強弱程度的物理量，能夠使實驗數據直接應用于寒區工程中，同時也可以和凍融指數進行對比研究，從而更進一步完善寒區研究。

因此，本文基于中國東北興安嶺地區18個氣象站點的地表溫度數據，對興安嶺地區年凍融頻次在空間上的分布情況和在時間上的變化情況進行調查研究，以期為之后在興安嶺地區的氣候變化、寒區工程、農林業生產發展等方面的研究提供可靠的數據支撐。

1 數據與方法

1.1 研究區概況

興安嶺地區有中國保存比較完好、面積最大的原始森林，森林覆蓋率達80.95%，區域范圍跨越黑龍江省和內蒙古自治區，屬于溫帶大陸季風性氣候。其中，大興安嶺平均海拔為573 m，冬季寒冷干燥且漫長，夏季溫暖濕潤而短暫，區域內河流較多，以黑龍江水系和嫩江水系為主。興安嶺地區是中國的北部邊境地區，包含極北部存在的不連續多年凍土區和一部分島狀多年凍土區以及南部的中-深季節凍土區，其凍土分布情況復雜，年凍融頻次變化差異較大，地理環境多樣。

本文選取興安嶺及其周邊部分地區作為研究區域。采用的氣象數據為興安嶺地區1990—2017年18 個氣象站0 cm 處逐日最高地表溫度、逐日最低地表溫度以及日平均地表溫度，圖1 為各站點分布情況，圖2為調查區凍土分布情況［16］。

圖1 興安嶺地區地形圖及氣象站點分布Fig. 1 Topographic map of Xing’anling region and distribution of meteorological stations

圖2 興安嶺地區凍土分布［16］Fig. 2 Distribution of seasonally frozen soil，island permafrost and discontinuous permafrost in Xing’anling region［16］

由于興安嶺地區各站點所處地理位置不同，氣候環境不同，各站點間年凍融頻次差異較大，全部站點在各年份的平均凍融頻次不能夠詳細地描述不同氣候環境下的年凍融頻次變化情況。研究區由南到北跨越中溫帶和寒溫帶兩個溫度帶，由西向東跨越半干旱區、半濕潤區和濕潤區三個氣候區。根據中國氣候區劃新方案（2010 年）［17］的東北地區氣候區劃結果對大興安嶺山脈沿線進行分區研究，如圖3 所示。寒溫帶濕潤區包含漠河、塔河、新林、圖里河等4個站點；中溫帶濕潤區包括呼瑪、加格達奇、孫吳、嫩江、黑河等5個站點；中溫帶半濕潤區包括額爾古納右旗、小二溝、博克圖、扎蘭屯、阿爾山、索倫等6個站點；中溫帶半干旱區包括滿洲里、海拉爾、新巴爾虎右旗等3個站點。

圖3 興安嶺地區氣候區劃［17］Fig. 3 Division of climate zone in Xing’anling region［17］

1.2 NDVI數據

本文采用的NDVI 數據為Land Term Data Record（LTDR）NDVI（AVH13C1）和MODIS NDVI（MOD13C2）兩種NDVI數據集的結合。其中1981—1999 年的NDVI 數據來自LTDR AVH13C1，該數據屬于日NDVI數據，其空間分辨率為0.05°，獲取數據的網址為http：//ltdr.nascom.nasa.gov/cgi-bin/ltdr/ltdrPage. cgi。2000—2014 年的NDVI 數據來自MODIS NDVI，MOD13C2屬于半月合成數據，其空間分辨率為0.05°，獲取數據的網址為https：//ladswe.nascom. nasa. gov。這兩種數據都已經經過幾何校正、大氣校正、輻射校正等預處理［18］，并且經過驗證，由LTDR 和MODIS 生成的1981—2014 年逐月NDVI數據可靠，可以用到本文的研究中［19］。

1.3 研究方法

1.3.1 年凍融頻次的定義

根據實驗室中對混凝土、巖石等寒區工程中涉及到的材料進行凍融循環實驗時所用到的凍融循環方法，為方便后續研究者們將實驗所得數據更好地應用到實際工程建設中，本文規定以一個測量日內日最高地表溫度大于0 ℃，同時日最低地表溫度小于0 ℃為一次凍融循環。一天內土壤由夜間最低溫度升至白天最高溫度時，經歷一次融化過程，再由白天最高溫度降至夜間最低溫度時，經歷一次凍結過程，視為一次完整的凍融循環。統計一個測量年內凍融循環的次數即為該年的年凍融頻次，以次/年表示。凍融循環次數由逐日觀測的日最高地表溫度和日最低地表溫度數據計算得出。

1.3.2 時空變化趨勢分析

在年凍融頻次時間變化趨勢的分析中，采用線性傾向估計的方法研究年凍融頻次在長時間序列上的變化趨勢及其傾向率。在年凍融頻次空間變化趨勢的分析中，基于ArcGIS 10.2，將多年平均凍融頻次和DEM 數據采用IDW（inverse distance weighted）插值方法進行插值，并采用線性擬合分析方法進一步研究年凍融頻次與經度、維度和海拔之間的變化關系。

1.3.3 突變分析

Pettitt 方法是通過檢驗序列均值變化的時間來確定序列突變點，是一種基于Mann-Whitney的非參數檢驗方法，可以對存在趨勢變化的序列進行突變檢驗。本文采用Pettitt檢驗法對興安嶺地區18個站點在長時間序列上的變化情況進行突變分析。在突變分析中需要用到的非參數統計量為

根據以上非參數統計量可以計算出

一般情況下，認為當P≤0.05時被檢驗序列中存在突變點。

1.3.4 權重分析

采用多元線性回歸權重分析的方法對興安嶺地區年凍融頻次在地理因素上的分布情況進行權重分析。在進行多元線性擬合之前，為去除經緯度、海拔數據在數值上對權重分析的影響，首先對經緯度和海拔數據進行標準化處理，然后對年凍融頻次與經緯度和海拔的標準化數據進行多元線性回歸擬合，整理出經緯度和海拔對年凍融頻次影響的權重。

2 結果與分析

2.1 年凍融頻次的時間變化

對興安嶺地區18 個站點的年凍融頻次數據求平均值，得到全部站點在各年份的平均凍融頻次。由圖4 可知，興安嶺地區年凍融頻次在1990—2017年間呈波動下降的趨勢，線性擬合曲線斜率為-2.0，即年凍融頻次平均每年下降2.0 次/年，通過了0.05 顯著性水平檢驗，年凍融頻次最高值達132 次/年，最低值為66次/年。

圖4 1990—2017年興安嶺地區年凍融頻次的年際變化Fig. 4 Interannual variation of annual freeze-thaw frequency in Xing’anling region from 1990 to 2017

圖5為4 個分區的年凍融頻次變化情況及其線性擬合線分布。從圖中可以看出，4 個分區內年凍融頻次的走勢與總體的年凍融頻次走勢相近，隨著時間的遷移，呈現明顯的下降趨勢，但各區年凍融頻次的高低和下降趨勢的大小不盡相同。圖5（a）表明，在興安嶺北部的寒溫帶濕潤區，年凍融頻次擬合曲線斜率為-2.0，即寒溫帶濕潤區年凍融頻次平均每年下降2.0 次/年；圖5（b）表明，在興安嶺東部的中溫帶濕潤區，年凍融頻次擬合曲線斜率為-1.8，即中溫帶濕潤區年凍融頻次平均每年下降1.8 次/年；圖5（c）表明，在興安嶺南部的中溫帶半濕潤區，年凍融頻次擬合曲線斜率為-2.3，即中溫帶半濕潤區年凍融頻次平均每年下降2.3 次/年；圖5（d）表明，在興安嶺西部的中溫帶半干旱區，年凍融頻次擬合曲線斜率為-1.9，即中溫帶半干旱區年凍融頻次平均每年下降1.9 次/年。以上線性擬合均通過了0.05顯著性水平檢驗。

圖5 興安嶺地區各氣候區年凍融頻次的年際變化Fig. 5 Interannual variation of annual freeze-thaw frequency of four climate zones in Xing’anling region：cold temperate humid area（a），middle temperate humid zone（b），middle temperate sub-humid zone（c）and middle temperate semi-arid zone（d）

注意到年凍融頻次在時間線上存在跳躍性較大的區域，因此為了進一步分析年凍融頻次在時間線變化中是否存在突變現象，本文采用Pettitt 法對興安嶺地區各站點在時間序列上的變化進行突變分析。對氣候區劃4個分區內各個站點的平均凍融頻次進行整合，得到各分區在每年的平均凍融頻次，之后對這一平均凍融頻次采用Pettitt 法進行整理，得到非參數統計量Ut，Nmax數據。由圖6 可知，興安嶺地區內4 個分區的統計量Ut，N的變化情況大致相同，統計量Ut，Nmax均在2004年達到最高值，從寒溫帶濕潤區到中溫帶半干旱區，各分區內的統計量Ut，Nmax分別為193、182、195、185，起對應的突變檢驗值P分別為1.08×10-4、3.20×10-4、8.77×10-5、2.39×10-4，所以各分區內檢驗值P均小于0.05，則可以說明興安嶺地區年凍融頻次在時間序列上的變化發生了顯著性突變，突變年份為2004年。

圖6 興安嶺地區各氣候區Ut，N的變化情況Fig. 6 Variation of Ut，N of four climate zones in Xing’anling region：cold temperate humid area（a），middle temperate humid zone（b），middle temperate sub-humid zone（c）and middle temperate semi-arid zone（d）

由于年凍融頻次在2004 年發生突變，所以以2004 年為分界線。1990—2004 年，興安嶺地區年凍融頻次平均在118 次/年左右，最大為132 次/年，最小為90 次/年；興安嶺北部的寒溫帶濕潤區年凍融頻次平均值在115 次/年左右，最大達到129 次/年，最小為96 次/年；興安嶺東部的中溫帶濕潤區年凍融頻次平均值在102 次/年左右，最大為114 次/年，最小為78 次/年；興安嶺南部的中溫帶半濕潤區年凍融頻次平均值在132 次/年左右，最大為152 次/年，最小為113 次/年；興安嶺西部中溫帶半干旱區年凍融頻次平均值在119 次/年左右，最大為151 次/年，最小為98次/年。

2005—2017 年，興安嶺地區整體多年平均凍融頻次在83 次/年左右，最大為95 次/年，最小為66 次/年；興安嶺北部的寒溫帶濕潤區年凍融頻次平均值在77 次/年左右，最大達到111 次/年，最小為59 次/年；興安嶺東部的中溫帶濕潤區年凍融頻次平均值在73 次/年左右，最大為97 次/年，最小為61 次/年；興安嶺南部的中溫帶半濕潤區年凍融頻次平均值在93 次/年左右，最高大為108 次/年，最小為75 次/年；興安嶺西部中溫帶半干旱區年凍融頻次平均值在88次/年左右，最大為102次/年，最小為68次/年。

2.2 年凍融頻次的空間變化

由于各站點年凍融頻次為長時間序列的數據，為了討論興安嶺地區各站點間年凍融頻次在空間上的變化情況，對各個站點的多年凍融頻次數據求平均值，得到各個站點1990—2017 年近30 年年凍融頻次的平均值。通過ArcGIS 10.2 對年凍融頻次數據進行處理，采用IDW 插值方法將年凍融頻次在空間上的變化情況展現在二維等值線地圖之上，與等高線插值圖進行對比研究，觀察年凍融頻次與海拔、經緯度之間的關系。圖7 為興安嶺地區平均凍融頻次插值圖，可以看出，興安嶺地區年凍融頻次整體上呈現西高東低、南高北低的趨勢。凍融頻次最大值出現在索倫站，平均凍融頻次為139 次/年；最小值出現在黑河站，平均凍融頻次為83次/年。

圖7 興安嶺地區多年平均凍融頻次的插值分布圖Fig. 7 Interpolation distribution of averaged annual freeze-thaw frequency in Xing’anling region

同時通過線性回歸分析的方法，對多年平均凍融頻次與海拔、經緯度進行線性擬合，更加精確地分析年凍融頻次在空間因素變化下的變化趨勢以及線性傾向率。如圖8（a）所示，對興安嶺地區各站點的平均凍融頻次與緯度進行線性擬合可以看出，興安嶺地區各站點間的平均凍融頻次在由南到北呈現不斷遞減的線性分布，確定系數R2為0.434，隨著緯度的增加，年凍融頻次大體上以5.5 次/度的速率遞減，擬合線通過了0.05 顯著性水平檢驗。圖8（b）為各站點的平均凍融頻次隨經度變化情況以及線性擬合，可以看出興安嶺地區各站點的平均凍融頻次由西向東呈現不斷遞減的線性分布，確定系數R2為0.284，隨著經度的增加，年凍融頻次大體上以2.4 次/度的速率遞減，擬合線通過了0.05 的顯著性水平檢驗。

如圖8（c）所示，興安嶺地區各站點的平均凍融頻次在垂直高度上隨海拔的升高呈現不斷遞增的線性分布，確定系數R2為0.140。隨著海拔的升高，年凍融頻次大體上以23 次/千米的速率遞增，擬合曲線未通過0.05 顯著性水平檢驗。斜率接近于0，說明年凍融頻次受海拔影響不大。

圖8 興安嶺地區多年平均凍融頻次隨經緯度及海拔的變化Fig. 8 Variation of averaged annual freeze-thaw frequency in Xing’anling region with latitude（a），longitude（b）and altitude（c）

由上述可知年凍融頻次在空間變化上會受到經緯度以及海拔的影響，并且隨緯度升高而降低，隨經度增加而降低，隨海拔升高而升高。為了進一步探究影響年凍融頻次在空間上變化的因素中經緯度和海拔的影響權重是多少，本文對經度、緯度和海拔在年凍融頻次的影響程度上利用多元線性回歸分析的方法進行權重分析。首先對經度、緯度和海拔的數據進行標準化處理，使原始數據去量綱化，從而消除原始數據在單位上的限制，使不同單位或量級的指標進行擬合，更好地進行權重分析。然后通過最小二乘法，以多年平均凍融次數為因變量，以經度、緯度和海拔的標準化數據為自變量，進行多元線性擬合，得到表1中的數據。

從表1 可以看出，在滿足0.05 顯著性水平檢驗的情況下，經度、緯度和海拔對興安嶺地區年凍融頻次多元線性擬合的擬合系數分別為-20.020、-28.580 和-6.332。如果僅考慮地理信息對興安嶺地區年凍融頻次的影響，對擬合系數進行標準化處理，得到經度、緯度和海拔對興安嶺地區年凍融頻次影響的權重分別為36.5%、52.0%和11.5%，由此可知在影響年凍融頻次空間變化的因素中，緯度影響最大，其次是經度，海拔影響最小。

表1 經緯度及海拔對興安嶺地區年凍融頻次影響的權重Table 1 Weight of influence of longitude，latitude and altitude on annual freeze-thaw frequency in Xing’anling region

2.3 年凍融頻次與NDVI變化情況對比

研究發現，北半球約70%的多年凍土地區在生長季（4—10 月）和72%的地區在秋季期間歸一化植被指數增加，凍土變暖與綠化呈正相關關系［20］。中國東北多年凍土區植被在生長季的NDVI 變化情況，在空間上大致表現為由西向東逐漸增加的趨勢，在時間上隨著時間的增加呈現顯著的增加趨勢。如圖9 所示，中國東北多年凍土區NDVI 由1990年的0.56增加到2014年的0.65。

圖9 1990—2014年東北地區NDVI的年際變化Fig. 9 Interannual variation of NDVI in Northeast China from 1990 to 2014

凍土區地表年凍融頻次除了會受到經緯度和海拔的影響之外，同時也受到植被覆蓋度、山坡坡向等多種因素的綜合影響。從圖10 可以看出，中國東北多年凍土區歸一化植被指數NDVI 隨時間序列呈顯著的波動上升趨勢，而年凍融頻次呈波動下降趨勢。采用Pettitt 法對東北地區NDVI 數據進行突變分析，得到NDVI 變化的統計量Ut，N在2000年取得最大值-150，其對應的突變檢驗值P=0.00049＜0.05，說明東北地區NDVI 在2000 年發生了突變。

圖10 興安嶺地區NDVI與年凍融頻次的年際變化對比Fig. 10 Comparison of interannual variation between NDVI and annual freeze-thaw frequency in Xing’anling region

圖11將大興安嶺山脈周邊地區的年凍融頻次與NDVI 變化情況進行相關性擬合，在擬合過程中令年凍融頻次滯后于對應NDVI 年份4 年，發現隨著植被覆蓋率的增加，興安嶺地區的年凍融頻次呈顯著下降趨勢，經線性擬合的確定系數達到0.460，由此可知興安嶺地區年凍融頻次與NDVI變化情況具有良好的負相關性。興安嶺地區年凍融頻次在2004 年發生突變，與NDVI 數據對比有明顯的滯后性，說明年凍融頻次變化會受到當地植被覆蓋率的影響，且這種影響有一定的滯后性。植被覆蓋率的增長會逐漸減緩地表對熱能的接收和散發，這個過程需要一定的時間，從而導致年凍融頻次的降低，這個過程需要一定的時間，所以導致年凍融頻次變化滯后于NDVI的變化。

圖11 興安嶺地區年凍融頻次與NDVI散點圖及線性擬合Fig. 11 Annual freeze-thaw frequency and NDVI scatter diagram and linear fitting in Xing’anling region

3 結論

本文通過對中國東北興安嶺地區及周邊地區18個站點1990—2017年的年凍融頻次進行整理，討論了興安嶺地區年凍融頻次在時間和空間上的分布情況，并利用權重分析、突變分析等方法，分析了地理因素對年凍融頻次變化的影響以及年凍融頻次的突變情況，并結合東北地區NDVI 變化情況進行了對比。得出以下結論：

（1）興安嶺地區年凍融頻次在空間上大體上呈現隨經度、緯度升高而降低，隨海拔升高而升高的趨勢，進一步分析得到，經度、緯度和海拔在影響年凍融頻次變化中以經緯度影響為主，海拔影響較小。

（2）受到氣候變化的影響，興安嶺地區年凍融頻次隨著時間的推移，整體呈現明顯下降的趨勢，且在2004 年前后發生突變。同時進行年凍融頻次與NDVI數據對比分析后發現，年凍融頻次與NDVI間有較好的負相關性，且東北地區NDVI 在2000 年左右發生突變，年凍融頻次的突變與之相比，存在一定滯后性，說明年凍融頻次在一定程度上受到當地植被覆蓋率變化情況的影響。

調研年凍融頻次分布情況可以使實驗室中得到的巖石凍融循環試驗數據直接對接寒區工程，對寒區工程的建設有著重要的指導意義，也可彌補凍融指數在工程實踐中的不足之處，同時對寒區的氣候變化和農林業發展等有著重要的借鑒意義。但本文數據存在空間分辨率和精度上的不足，在討論年凍融頻次的影響因素時，也只考慮了經緯度和海拔的影響，對其他影響因素暫未考慮，期待未來加深研究。