寒凍風化控制的祁連山風化碎屑的空間分布

洪 洋， 耿豪鵬，2， 潘保田，2

（1.蘭州大學資源環境學院西部環境教育部重點實驗室，甘肅蘭州 730000；2.甘肅省石羊河流域野外科學觀測研究站，甘肅蘭州 730000）

0 引言

地貌是由巖體的抬升、變形和風化，以及地表的侵蝕、搬運和堆積等過程所塑造的［1］，其中風化作用是地貌侵蝕過程的先導［2］。寒凍風化是冰緣地貌區及高寒山區的主要物理風化過程［3-5］。20世紀80年代起，崔之久等通過野外定位觀測探討了寒凍風化的主要特點及影響因素［6-7］，發現寒凍風化提供的風化碎屑物質是冰緣地貌區及高寒山區崩塌、泥石流等災害性地貌過程［8］及倒石堆的重要物質來源［9］，寒凍風化可以通過影響倒石堆的運動速率，進而破壞公路等基礎設施［10］，由此開啟了我國研究寒凍風化過程的序幕。

寒凍風化的作用機制存在兩種不同的假說。第一種是凍融作用（freeze-thaw）相關的體積膨脹理論，即巖石孔隙中的水從液態變為固態時體積膨脹，對巖壁施加壓力使巖體破碎，易于崩解［11］。凍融作用要求存在一個飽和且密閉的空間，使體積膨脹所產生的壓力可以全部施加于巖壁［12］。但自然界中基本不存在飽和且密閉的孔隙，不飽和孔隙中體積膨脹所產生的壓力僅為飽和孔隙的10%左右［13］，因此凍融作用對真實地貌的有效性一直存在爭議［14］。為此學者提出分凝冰理論（segregation ice growth）來解釋真實環境中的寒凍風化機制。分凝冰理論認為高寒環境下巖石孔隙中廣泛存在分凝冰，由于界面融化和曲率效應，過冷水在冰-巖界面間形成薄膜［15］，使液態水從溫度較高的區域流向分凝冰所在孔隙，促進分凝冰不斷生長［16］，冰-巖界面間的范德華力及靜電力使巖體破碎［17］。分凝冰理論在實驗室得到了廣泛的驗證［18-20］，以分凝冰理論為基礎構建的寒凍風化模型也廣泛應用于計算不同區域的寒凍風化強度，進而解釋高寒山區侵蝕速率的時空變化及現代地貌單元的空間分布特征［21-23］。因此本文中的寒凍風化指基于分凝冰理論的巖體破碎過程。

廣泛分布的風化碎屑是高寒山區及冰緣地貌區寒凍風化的主要產物［24］，極易受多種因素誘發轉化為滑坡、泥石流等地質災害［25］。因此探明風化碎屑空間分布的主要控制因素對預防地質災害的發生具有重要意義。由于分凝冰生長速率及液態水從高溫向低溫的流動過程均與溫度相關［12，17］，因此溫度是寒凍風化的主要控制因素。同時寒凍風化對巖體破碎的影響是緩慢而持續的累積過程，因此在探討寒凍風化對現代地貌單元空間分布特征的影響及其在長期地貌演化中的作用時考慮過去溫度變化（包含一個或多個冰期—間冰期旋回）的影響十分必要［26-28］。同時前人研究表明，末次冰盛期以來的溫度變化計算的寒凍風化強度可以很好地約束風化碎屑空間分布的下界［29-30］。但上述研究均在新西蘭南阿爾卑斯山、意大利阿爾卑斯山等較為濕潤的山區開展，干旱半干旱高寒山區溫度變化時間尺度的選取尚需進一步探究。

前期野外考察結合遙感影像發現祁連山北部風化碎屑物質廣泛分布，且存在由風化碎屑轉化產生的滑坡、泥石流等地質災害。因此本文擬以祁連山北部為研究區，探討干旱半干旱高寒山區風化碎屑空間分布的主要控制因素。基于以上目標，我們首先利用哨兵二號遙感影像，采取NDVI（歸一化植被指數）閾值、海拔篩選結合目視解譯的方法提取祁連山北部風化碎屑的空間分布；其次利用寒凍風化模型計算不同時間尺度溫度變化情景下的寒凍風化強度的空間分布特征；最后對比分析風化碎屑及寒凍風化強度的空間相關性，探討干旱半干旱高寒山區風化碎屑空間分布的主要控制因素。

1 研究區概況

祁連山位于青藏高原東北緣，晚新生代起開始快速抬升［31］。本研究所關注的祁連山北部海拔較高（平均海拔約3 700 m），坡度較陡（超三分之一區域坡度超過20°），年平均氣溫約2.3℃，年降水量在100～600 mm·a-1之間變化，且呈現明顯的東高西低的分布趨勢［32］。植被表現出明顯的垂直地帶性：2 000～2 300 m主要為荒漠草原帶，2 300～2 600 m主要為溫性草原帶，2 600～3 200 m主要為森林草原帶，3 200～3 700 m主要為灌叢草原帶，3 700～4 000 m主要為高寒草甸帶，4 000 m以上主要為冰雪帶［33］；同時植被還表現出明顯的水平差異，植被垂直帶譜由東向西趨于簡化［34］。祁連山北部是疏勒河、黑河及石羊河三大內陸河流域的發源地，其中共有58條河流的多年平均出山徑流量大于0.1 m3·s-1［32］。本文選取祁連山北部烏鞘嶺以西為研究區（圖1）。

圖1 研究區數字高程模型（Shuttle Radar Topography Mission DEM，http：//srtm.csi.cgiar.org/）及氣象站點位置（黃色圓點）Fig.1 The digital elevation model（Shuttle Radar Topography Mission DEM，http：//srtm.csi.cgiar.org/）of the study area with locations of meteorological stations（yellow circles）

2 數據和方法

2.1 風化碎屑邊界的提取

我們采取NDVI（歸一化植被指數）閾值、海拔篩選結合遙感影像目視解譯的方式進行風化碎屑邊界的提取。風化碎屑與植被由于反射率不同在遙感影像上表現出的光譜特征差異是我們采用NDVI閾值提取與目視解譯的基礎。我們使用哨兵二號10 m分辨率的遙感影像（European Space Agency，https：//sentinel.esa.int/web/sentinel/home）用于計算研究區的NDVI，所用影像選取原則為2021年6—9月以最大限度去除積雪的影響，且云覆蓋量低于10%。由于植被垂直帶譜在東西方向存在差異，因此NDVI閾值的選取隨植被特征的變化而變化，范圍為0.09～0.15。谷歌地球0.6 m分辨率的遙感影像主要用于目視解譯，對NDVI閾值提取結果進行修正，去除水體、云霧等因素干擾。同時為了剔除低海拔區域其他成因裸巖斑塊的影響，我們依據研究區植被的垂直地帶性特征結合野外考察、遙感影像目視解譯及前人研究成果［35］，以風化碎屑萎縮或擴張變化的海拔下限3 800 m為閾值對提取結果進行篩選。我們利用30 m分辨率的數字高程模型（Shuttle Radar Topography Mission DEM，http：//srtm.csi.cgiar.org/）進行研究區海拔屬性的提取。

2.2 寒凍風化模型

Hales等［36］以分凝冰理論為基礎設計了寒凍風化模型，以寒凍風化強度（frost cracking intensity，Fci）作為評價溫度對巖體破壞程度的指標。依據Hales等［36］的模型，我們先用一維熱傳導公式模擬溫度隨時間（t）和深度（z）的變化（圖2）：

式中：T是不同時間（t，單位為天）不同深度（z，單位為cm）下的溫度（℃）；TMA是年平均氣溫（℃）；Ta是年氣溫變幅的一半（℃）；α是熱擴散系數（設定為1 mm2·s-1）；Py指氣溫變化的周期，在本研究中為1年。寒凍風化模型的構建存在三個假設前提：（1）當溫度在-3～-8℃時（即寒凍風化窗）分凝冰生長速率最快［12］；（2）當溫度高于0℃時存在液態水；（3）液態水從高溫向低溫單向流動。依據Hales等［36］的定義，寒凍風化強度（Fci，°C·cm-1）是寒凍風化窗內溫度梯度隨時間的積分［37］：

式中：?T是溫度梯度；n是符合三個假設前提的天數。同時，由于研究區屬于干旱半干旱區，還需考慮液態水在流動過程中因無法及時補充造成的可能損耗，因此我們引入Anderson等［27］提出的距離函數（e-D/D*），以Fci與距離函數的乘積表征適用于干旱半干旱區的寒凍風化強度［37］。通過對修改后的寒凍風化強度進行深度積分，得到不同年平均氣溫條件下的總寒凍風化強度。依據氣溫和海拔之間的線性關系可得到不同年平均氣溫所對應的海拔，即可得到不同海拔條件下的總寒凍風化強度［37］：

式中：D是0℃所在位置到-3℃所在位置的距離；D*是長度指數，依據Anderson等［27］設置為50 cm。但研究發現某些日期溫度隨深度的變化并非單調的，而是先增加后減小（圖2紫色實線），因此存在一個溫度轉折點，導致0℃同時出現在兩個不同的深度。因此在本研究中，將D設定為轉折點所在位置到-3℃所在位置的距離以保證溫度隨深度的單調變化。

圖2 年平均氣溫為-3℃及Ta=15℃時溫度隨時間和深度的變化［灰色陰影為寒凍風化窗所在范圍（溫度在-3～-8℃時），垂直的黑色虛線為0℃所在位置；紫色實線（t=182）代表溫度隨深度的非單調變化］Fig.2 Schematic representation of the temperature variation as a function of depth using TMA of-3°C and Ta of 15°C［The grey shadow represents the frost cracking window，between-3 and-8°C；The black dashed line represents the location of 0°C；The purple line（t=182）is an example of a depth profile with two 0°C locations and an inflection point］

為了計算平均寒凍風化強度來評估不同時間尺度溫度變化對寒凍風化的累積影響，我們從Vostok冰芯記錄中以500年為時間間隔提取不同年代的年平均氣溫（TMA）并假定Ta不變，得到不同年代的總寒凍風化強度，取其平均值即為平均寒凍風化強度［37］：

式中：i是時間間隔的數量；Fci-ti是第i個時間間隔相對應年代的總寒凍風化強度。

2.3 溫度數據

年平均氣溫可以通過11個氣象站點獲得（圖1中的黃色圓點，表1），數據來源為國家氣象信息中心（https：//data.cma.cn/），記錄了1981—2010年的日平均氣溫的平均值。我們利用氣象站記錄的氣溫與站點海拔之間的線性回歸關系來計算氣溫直減率及年平均氣溫0℃所在的海拔高度（圖3）。我們利用中國區域地面氣象要素驅動數據集（1979—2018）空間分辨率為0.1°，時間分辨率為3小時的氣溫數據［38-39］計算年氣溫變幅（日均最高氣溫與最低氣溫的差值，其數值的一半即Ta，圖4）。由于古里雅冰芯及敦德冰芯的年代存在爭議［40］，南極Vostok冰芯［41］所記錄的末次冰盛期（Last Glacial Maximum，LGM）的降溫幅度與青藏高原6～9℃的降溫幅度較為一致［42］，其用于計算LGM以來平均寒凍風化強度的有效性在祁連山東段得到了驗證［37］，且其溫度記錄的時間范圍為420 ka至今，滿足本研究對不同時間尺度溫度變化的需求，因此我們選擇南極Vostok冰芯作為長期溫度變化的參考（圖5）。

表1 11個國家氣象信息中心氣象站基本信息表Table 1 Information of 11 national meteorological stations

圖3 11個氣象站點海拔與年平均氣溫散點圖Fig.3 Scatter plots of elevation against TMA

圖4 祁連山Ta空間分布圖Fig.4 The spatial distribution of Ta in the Qilian Mountains

3 結果與分析

3.1 風化碎屑邊界的空間分布特征

風化碎屑的范圍表現為自東向西逐漸增大的趨勢（圖6）。為了進一步驗證風化碎屑邊界提取的準確性，我們在祁連山東中西段各選擇一個典型區域對比提取的風化碎屑邊界與哨兵二號遙感影像［圖6（b）～6（d）］。結果顯示祁連山東段［圖6（b）］與中段［圖6（c）］的風化碎屑邊界與遙感影像所顯示的風化碎屑邊界具有較好的一致性。而祁連山西段年降水量較低，植被覆蓋較差（NDVI均值0.05），且多為斑塊狀分布，導致利用NDVI閾值提取風化碎屑邊界較為困難，因此風化碎屑邊界與遙感影像所顯示的風化碎屑邊界存在部分區域的不匹配［圖6（d）］。

圖6 祁連山風化碎屑邊界分布圖［（b）～（d）是黑色實線框區域內哨兵二號遙感影像及風化碎屑邊界］Fig.6 The spatial distribution of the boundary of scree in the Qilian Mountains：（b）～（d）are Sentinel-2 imageries and the boundary of scree in the range of black rectangle

3.2 寒凍風化強度的時空變化特征

25～16 ka（圖7藍色實線）及16 ka以來（圖7紫色實線）平均寒凍風化強度隨海拔的變化曲線表現為單峰特征。由于從16 ka開始溫度快速升高且在一段時間內維持在相對高溫狀態（圖5虛線框），寒凍風化窗向高海拔區域移動，導致16 ka以來平均寒凍風化強度的峰值海拔高于25～16 ka，且兩者寒凍風化峰值海拔差距較大，因此25 ka以來（包含一個冰期-間冰期旋回）平均寒凍風化強度隨海拔的變化曲線表現為雙峰特征（圖7橙色實線）。十萬年尺度平均寒凍風化強度隨海拔的變化表現為單峰特征（圖7），最大值在8～9℃·cm-1之間，且均在海拔3 900～4 700 m之間達到最大值（圖7灰色陰影），不受隨機選取的時間范圍的影響。其與25 ka平均寒凍風化強度隨海拔變化的雙峰特征存在較大差異，主要表現為海拔5 000 m以上，25 ka平均寒凍風化強度依然較高。而在海拔5 000 m以下，25 ka平均寒凍風化強度隨海拔的變化趨勢與十萬年尺度平均寒凍風化較為一致。

圖5 420 ka以來Vostok冰芯溫度記錄（修改自Petit等［41］，紫色虛線框為25 ka以來溫度變化）Fig.5 420 ka continuous temperature record based on the Vostok ice core，Antarctica：modified from Petit et al［41］.Purple dashed box is 25 ka continuous temperature record

圖7 不同時間尺度平均寒凍風化強度隨海拔的變化（灰色陰影代表十萬年尺度平均寒凍風化強度峰值的海拔范圍）Fig.7 The variation of time averaged Fci with elevation for different temporal scale（The grey shadows represent the peak elevation range of time averaged frost cracking intensity）

3.3 年氣溫變幅對寒凍風化強度空間變化的影響

為了探討Ta對寒凍風化強度空間變化的影響，我們計算了不同Ta條件下總寒凍風化強度隨海拔的變化［圖8（a）］，并以25 ka平均寒凍風化強度為例，計算了不同Ta條件下平均寒凍風化強度隨海拔的變化［圖8（b）］。結果顯示，Ta的增大可以使相同海拔條件下總寒凍風化強度增加，同時總寒凍風化強度峰值所在的海拔范圍增大，進而導致平均寒凍風化強度增加，同時其峰值海拔范圍增大。因此本文采取的差異化Ta參數是準確計算研究區寒凍風化強度的關鍵。

圖8 不同Ta條件下總寒凍風化強度（a）及25 ka平均寒凍風化強度（b）隨海拔的變化（灰色陰影代表Ta=15℃時寒凍風化強度峰值的海拔范圍）Fig.8 The variation of total Fci（a）and time averaged Fci（b）with elevation with different Ta（The grey shadows represent the peak elevation range of Fci with Ta=15℃）

3.4 風化碎屑邊界與寒凍風化強度的空間關系

由于萬年尺度與十萬年尺度平均寒凍風化強度的空間分布特征存在差異（圖7），因此我們分別制作了風化碎屑邊界與萬年尺度［25 ka，圖9（a）］和十萬年尺度［100 ka，圖9（b）］平均寒凍風化強度的空間分布圖，并以風化碎屑與寒凍風化強度高值重疊面積占風化碎屑面積的百分比（表2）來定量探討風化碎屑邊界與不同時間尺度溫度變化控制的平均寒凍風化強度的空間關系。結果表明，25 ka和100 ka平均寒凍風化強度的空間分布與風化碎屑邊界的空間分布均存在較好的相關性，重疊面積占比大于65%，其中祁連山西段的重疊面積可達90%以上。由于Ta呈自東向西增大的趨勢（圖4），平均寒凍風化強度的峰值海拔范圍增大［圖8（b）］，導致祁連山西段寒凍風化強度高值的覆蓋面積占比（大于65%）大于祁連山東段（小于50%）。但祁連山西段100 ka平均寒凍風化強度高值覆蓋面積占比（76.49%）較25 ka平均寒凍風化強度高值覆蓋面積占比（65.89%）的顯著增加（超過10%）并沒有導致重疊面積占比的大幅增加（僅增加3.4%），表明100 ka平均寒凍風化強度高值面積的增大并沒有很好地提高其對風化碎屑邊界的約束，而25 ka平均寒凍風化強度的高值以更小的覆蓋面積約束了90%的風化碎屑范圍，即25 ka平均寒凍風化強度的高值與風化碎屑邊界的空間相關性優于100 ka。可能原因為祁連山西段海拔5 000 m以上區域所占比重較大（約1.5%），而25 ka平均寒凍風化強度在海拔5 000 m以上依然較高（圖7橙色實線），可以較好地約束風化碎屑的邊界。而對于祁連山中東段，100 ka平均寒凍風化強度高值覆蓋面積占比（47.21%）較25 ka平均寒凍風化強度高值覆蓋面積占比（27.78%）的顯著增加（約20%）導致重疊面積占比的大幅增加（約20%），表明25 ka和100 ka平均寒凍風化強度的高值均可以很好地約束風化碎屑邊界。可能原因為祁連山中段及東段海拔較低（99.8%的區域海拔低于5 000 m，圖1），在海拔5 000 m以下，25 ka和100 ka平均寒凍風化強度隨海拔的變化趨勢較為一致（圖7）。

表2 風化碎屑及不同時間尺度平均寒凍風化強度面積統計Table 2 Statistics of area of scree and different temporal of time averaged Fci

圖9 祁連山25 ka（a）及100 ka（b）平均寒凍風化強度及風化碎屑邊界空間分布圖Fig.9 The spatial distribution of time averaged Fci since 25 ka（a）and 100 ka（b）and the boundary of scree in the Qilian Mountains

4 討論

萬年尺度平均寒凍風化強度隨海拔的變化趨勢在5 000 m以下與十萬年尺度平均寒凍風化一致，表明寒凍風化在萬年及以上時間尺度對該區域的作用是持續且穩定的。這種持續且穩定的寒凍風化為滑坡、泥石流等災害性地貌過程提供充足的風化碎屑物質。因此在探討寒凍風化對長時間（萬年及以上）尺度地貌演化的累積影響時，時間尺度的選取對結果的影響可以忽略不計。我們推測的可能原因為：在冰期-間冰期旋回中，冰期持續時間長，間冰期持續時間短（圖5），因此當時間尺度足夠長時（例如萬年及以上），相對短暫的間冰期在累積寒凍風化中的影響相對于冰期可以忽略不計。因此，在探討寒凍風化在長期地貌演化中的累積影響時，如果研究區溫度記錄無法覆蓋研究所需的較長時間范圍（如百萬年），可以采用萬年或十萬年尺度溫度記錄計算的平均寒凍風化強度作為替代。

25 ka平均寒凍風化強度與祁連山西段風化碎屑邊界的空間相關性優于100 ka平均寒凍風化強度，表明當探討寒凍風化對現代地貌單元空間分布特征的影響時，間冰期的快速升溫過程顯得尤為重要。溫度快速升高且在一段時間內（千年—萬年）維持在相對高溫狀態（圖5虛線框），導致寒凍風化窗向高海拔區域移動，形成萬年尺度平均寒凍風化強度隨海拔變化曲線中的高海拔峰值（圖7）。但由于快速升溫為非連續性過程，寒凍風化對該峰值海拔范圍的影響是階段性的。在冰期大規模冰進的背景下，前一個間冰期所積累的碎屑物質無法保存，即寒凍風化的累積影響近乎清零。因此在探討寒凍風化對現代地貌單元空間分布特征的影響時應采用最近一次冰期-間冰期旋回的溫度變化數據（即末次冰盛期）。

在全球變暖的背景下，前人對由凍土退化、冰川退縮、冰湖擴張等引發的災害性過程的研究已較為成熟［44-46］，但對受寒凍風化影響的災害性地貌過程缺乏重視。本研究表明，溫度升高使寒凍風化的作用區向高海拔區域移動，其提供的風化碎屑物質是滑坡、泥石流等災害性地貌過程重要物質來源（圖10），進而產生新的高風險災害區域。寒凍風化不僅可以通過控制滑坡、泥石流等災害性地貌過程的發生，成為滑坡-堰塞-潰決洪水這一災害鏈［47］的源頭，還可以加速與冰川退縮相關的災害性地貌過程的發生［48］。本研究強調了全球變暖背景下寒凍風化對災害性地貌過程的調控，為全球變暖背景下災害性地貌過程的預測提供了新的思路，可以成為防災減災決策的重要參考。本研究是對崔之久先生寒凍風化如何影響生產建設并服務于社會經濟建設的繼承與發展。

圖10 祁連山典型風化碎屑積累實體例證Fig.10 The accumulation of scree in the Qilian Mountains

5 結論

本文通過NDVI閾值、海拔篩選結合遙感影像目視解譯，提取了祁連山北部的風化碎屑邊界，通過對比分析風化碎屑邊界及不同時間尺度平均寒凍風化強度的空間相關性，探討干旱半干旱高寒山區風化碎屑空間分布的主要控制因素。獲得主要結論如下：

（1）萬年尺度平均寒凍風化強度隨海拔的變化曲線表現為雙峰特征，十萬年尺度平均寒凍風化強度隨海拔的變化曲線表現為單峰特征。海拔5 000 m以下萬年尺度與十萬年尺度平均寒凍風化強度隨海拔的變化趨勢一致。

（2）在祁連山北部萬年尺度及十萬年平均寒凍風化強度的高值與風化碎屑邊界均存在較好的空間一致性，祁連山西段萬年尺度平均寒凍風化強度的高值與風化碎屑邊界的空間相關性優于十萬年尺度。

（3）探討寒凍風化對長時間（萬年及以上）尺度地貌演化的累積影響時，時間尺度的選取對結果的影響可以忽略不計；探討寒凍風化對現代地貌單元空間分布特征的累積影響時，應采用最近一次冰期-間冰期旋回的溫度變化數據。

（4）在全球變暖背景下，應重視寒凍風化對災害性地貌過程的調控。