高原寒區地表土壤凍融作用的空間參數化表征研究

高會然， 許 沖， 張萬昌， 易亞寧， 肖子亢

（1. 應急管理部 國家自然災害防治研究院，北京 100085； 2. 中國科學院 空天信息創新研究院，北京 100094；3. 復合鏈生自然災害動力學應急管理部重點實驗室，北京 100085）

0 引言

在全球氣候持續變暖的背景下，冰凍圈領域已成為研究熱點［1-4］。當前，冰凍圈研究不僅關注其在氣候系統中的作用，還更加關注對人類社會可持續發展的現實影響和潛在威脅［5-7］。作為冰凍圈的重要組成部分，凍土的變化或退化對生態環境系統、人類生活和生產安全等方面的直接或間接風險，這些研究方向也越來越受到學者的關注［8-9］。

當前，凍土在全球范圍內正呈區域性退化趨勢，主要表現為多年凍土面積縮小、活動層加厚，季節凍土層變薄、凍結時長縮短等［10-13］。與多年凍土變化相比，季節凍土凍融循環過程周期短、范圍廣，對陸地表面與大氣之間的物質與能量交換、陸地表面景觀格局的影響劇烈而深遠［5-6，14-16］?，F有研究表明，季節凍土對地表景觀的影響機制主要體現在凍土凍融的風化作用、凍脹和融沉作用以及凍結滯水作用［17］。一方面，季節凍土是深度不斷變化的隔水層，并通過影響區域土壤水特性而形成特殊的水文過程和生態過程［18-19］。另一方面，在地表季節凍土凍融過程中，山體邊坡土壤結構、坡體的穩定性、空隙水壓力、土體摩擦力等物理特性發生不同程度的變化，為淺層地質災害的發育提供致災條件和驅動力［20-23］。因此，摸清地表季節凍土的分布格局及其時空變化對開展寒區自然科學研究、保障或維護人類的生產活動和生存環境安全均具有重要的研究意義。

傳統的凍土研究多采用地面觀測的方法，但是傳統方法在數據實時獲取、信息快速提取以及特征參量空間離散化的研究需求下面臨重大挑戰［24-26］。20世紀80年代以來，國內外學者利用遙感技術進行了大量的區域或全球尺度的凍土時空動態遙感監測研究，并取得一系列數據或技術成果［27-32］。但是，遙感凍土監測受限于當前遙感技術水平，低空間分辨率的遙感凍土數據產品在中小區域尺度上的實際應用仍十分有限。隨著凍土物理學的發展以及對凍土系統中的水熱傳輸特性和機制的深入研究，國內外學者發展出了多種基于物理機理的凍土水熱傳輸過程數值模型。例如，Flerchinger 等［32］構建的一維水熱耦合過程模型SHAW（Simultaneous Heat and Water）模型，Wood等［33］提出的適用于大區域尺度水文過程模擬的VIC（Variable Infiltration Capacity）模型，陳仁升等［34］建立的嵌入分布式水文過程模型的高寒山區水熱耦合模型（DWHC），Gao等［31］基于水熱耦合原理研發了一套空間全分布式凍土水熱過程數值模型（Fully Distributed Frozen Soil Hydrothermal Processes Integrated Modeling System，FFIMS 模型），以及近十年來提出或改進的多種考慮凍融循環的陸面過程模型，如GEOTOP 模型、CLM 模型、UW-VIC 模型等［35-38］，為凍土系統動態過程研究提供了有效的方法和手段。

當前，基于物理機理的凍土水熱傳輸過程數值模型研究取得了較大的進展，部分陸面過程模型通過模型集成的方式，具備了精細刻畫地表土壤凍融循環多參量和多過程的能力［39-40］。但是，目前凍土變化監測與模擬研究成果在陸面環境要素響應和凍土災害防治等領域的應用能力仍處于較低的水平。一方面，高時空分辨率凍土特征參量數據獲取難度相對較大，限制了大范圍長時段凍土凍融循環及其凍融作用的定量分析；另一方面，地表土壤凍融作用作為凍土水熱過程與陸面景觀格局演變相互聯系的重要紐帶，其表征指標識別與空間參數化的定量研究不足。區域尺度凍土凍融作用的參數化表征及其進一步在災害風險評估與防控領域中的應用研究仍是空白。

因此，為了精細刻畫全球氣候變化下的季節凍土水熱過程，系統性表達凍土凍融循環的時空變化特征以及凍融作用特征，本研究以地處青藏高原東南緣的高山峽谷區及其周邊地區為研究區，綜合考慮氣象、植被、積雪與土壤等多個陸面過程以及各個過程中的水熱動態平衡，在基于凍土水熱耦合系統性理論的空間全分布式的凍土水熱耦合過程數值模型FFIMS 模型框架下，進行了適用于青藏高原高海拔凍土區的凍土水熱過程數值建模。在研究區開展長時段（2010年8月1日至2020年7月31日，共10 個水文年）逐日凍土水熱過程模擬，獲取研究時段內高精度的凍土系統空間分布式水熱過程參量，并著重分析研究區凍土環境特征及其時空變化特征。在此基礎上，結合凍土邊坡穩定性理論，創新性提出凍土凍融作用的空間參數化方法，定量研究地表土壤凍融循環對土體結構穩定性的影響，為區域尺度凍土相關研究和凍融災害風險評估與災害防控等研究提供思路借鑒和數據及技術支撐。

1 理論與方法

1.1 總體方法框架

空間分布式的季節凍土特征參數是開展地表土壤凍融作用空間參數化的數據基礎，本研究首先利用現有的基于物理機制的分布式凍土水熱過程數值模型，獲取研究區長時間序列的凍土特征參數數據集，包括植被冠層、積雪、土壤剖面等凍土系統多要素、多過程空間水熱參量。然后，利用GIS空間分析、數據挖掘等數據處理技術，分析近十年青藏高原東南緣季節凍土凍融循環及其變化特征，揭示氣候變化下的復雜地理環境中的凍土時空演化規律。最后，根據土壤剖面水熱過程的多種凍土特征參量，結合土體抗剪強度和凍土邊坡穩定性相關理論，提出可以表征凍土凍融作用的空間參數化方案。主要研究框架如圖1所示。

圖1 本研究的主要方法框架示意圖Fig. 1 Overall framework of the integrated approach adopted in this study

1.2 分布式凍土水熱過程數值建模

作為凍土過程數值模型基礎理論的水熱耦合原理考慮了氣象、植被冠層、積雪與土壤等多個過程以及各個過程中的水熱平衡［31，34，41-42］，各個過程中的水熱平衡方程見式（1）～（5）。

（1） 植被冠層

植被及其冠層設計為單層單節點結構，同時考慮了植被生長和季節變化特征，冠層能量與水平衡方程如式（1）和式（2）所示。

式中：ρa、ca和T分別為大氣密度（kg·m-3）、大氣比熱容（J·kg-1·℃-1）和冠層溫度（℃）；ke為冠層熱量傳導系數（m2·s-1）；Hl為冠層向大氣傳輸的熱量（W·m-3）；ρv為蒸氣密度（kg·m-3）；El為冠層蒸散發量（mm）。

（2） 積雪

將積雪設計為單層雙節點的結構，考慮積雪的累積和升華以及水熱在雪層中的傳導過程，積雪能量平衡方程可用式（3）表示。

式中：ρsp、wsp和ksp分別為雪密度（kg·m-3）、體積含水量（m3·m-3）和熱傳導系數（W·m-1·℃-1）；ci和ρl分別為冰的比熱容（J·kg-1·℃-1）和水密度（kg·m-3）；Rn為積雪下的凈輻射通量（W·m-2）；Lf和Ls分別為融化潛熱與升華潛熱（J·kg-1）；qv為水汽通量（kg·s-1·m-2）。

（3） 土壤

設計為不同網格不同深度的多層土壤結構，考慮了凍融過程對水熱通量的影響，土壤能量和水平衡方程可分別表示為式（4）和式（5）。

式中：Cs為土壤比熱容（J·kg-1·℃-1）；ks為土壤熱傳導系數（W·m-1·℃-1）；K為土壤導水率（m·s-1）；ψ為土壤水勢（m）；U為土壤水通量的源/匯項（m3·m-3·s-1）。

在FFIMS 模型的水熱平衡方程的求解過程中，將植被冠層、積雪和土壤剖面分為有限層數，每層用一個節點表示，每個節點中的各水熱分量的存儲量以各節點所在層的土壤厚度為準。利用隱式有限差分方程組求解每個節點以及節點間的能量和水的平衡方程［43］。模型中各個子過程的理論和方法涉及大量參數和變量，包括一些可調參數、經驗常數和物理常數等。其中，可調參數主要是指受研究區域環境條件影響的參數，需要根據研究區地形、氣象、生態等條件調整。本研究建立的模型中主要的參數取值及其說明見表1。

表1 FFIMS模型可調參數及其在本研究中的取值Table 1 Variable parameters of FFIMS model and their values in this study

FFIMS 模型的驅動數據、輸入和輸出均為空間分布式的數據，本研究中所有的數據均統一采樣為1 km 空間分辨率，模型模擬的時間步長為24 h（逐日模擬）。其中，模型驅動數據和輸入數據主要包括氣象驅動數據（氣溫、降水、風速、相對濕度、大氣壓強、日照時長等要素）、土壤屬性數據（土壤類型、土壤深度、土壤機械組成、有機質含量等屬性）、地表覆被屬性數據（葉面積指數、葉面特征、植被高度、干生物量等）、土壤溫濕度初始狀態數據。模型輸出數據主要是可以表征凍土系統特征或狀態的變量，包括冠層水熱參量（冠層輻射、冠層截留等）、積雪水熱參量（雪面輻射、積雪厚度、雪水當量等）、土壤剖面水熱參量（包括土壤剖面溫度、土壤剖面含水/冰量、地表熱通量等）。

1.3 凍融作用參數化表征方法

地表土壤的凍融作用直接影響或引起土體分裂以及土壤顆粒的重新組合，從而影響土體的抗剪強度，導致土體從穩定狀態轉向相對不穩定狀態，威脅地區的邊坡穩定、工程安全等［44-45］。但是，凍融作用對土抗剪強度以及相關指標的影響程度目前仍無定論。根據本研究獲取的精細凍土特征參量，結合現有的理論研究基礎，可以實現定量化表征凍土凍融作用對土體抗剪強度指標的相對損傷程度。因此，本研究以土壤剖面含水量和土壤剖面含冰量兩個關鍵凍土特征參量為主要參數，構建土體抗剪強度的損傷系數。

土壤內摩擦角和黏聚力是表征土壤抗剪強度的兩個關鍵指標。已有研究表明，土體抗剪強度隨著土壤含水率的增加而降低，其中土壤內摩擦角隨土壤含水率增加而減小，黏聚力隨之多呈現先增大后減小的趨勢［46-47］。土壤在凍結狀態下，土體中既有固體冰也有液態水，此時的土體抗剪強度主要由土體固有特性和固體冰以及兩者間的相互作用決定的［48-50］。一般情況下，土壤含冰量越高，土體抗剪強度越大。因此，本研究以統計時段內的土壤剖面含冰量的變化量為依據，將土體抗剪強度損傷系數表示為e的指數函數，見式（6）。

式中：f為土體抗剪強度損傷系數（0≤f＜1），f值越小，表示凍融作用對土體抗剪強度的損傷越大；Δθice為某時段內土壤含冰量的變化量（m3·m-3），計算方法見式（7）；k為關于土壤含水量的系數，計算方法見式（8）。

式中：D為統計時間段（d）；θice，d為d時刻的土壤剖面含冰量（m3·m-3）；θw，d為d時刻的土壤剖面含水量（m3·m-3）；θ'D為D時段土壤剖面含水量的平均值（m3·m-3）。

式（6）～（8）表示，隨著凍土融化，土壤含冰量降低，土壤含水量上升，土壤含冰量的變化量增加，f值隨之下降，意味著凍融作用對土體穩定性的損傷加大。土壤含水量通過影響土壤黏聚力從而對土體抗剪強度產生影響，關于土壤含水量的調節系數k決定了f值變化的程度。當土壤開始融化至融化初期，在土壤含水量低于θ'D時，土壤含水量緩慢增加，k值不斷增大，f值下降緩慢甚至有所增加，此時土體抗剪強度由土壤水和固體冰以及兩者間的相互作用等條件決定；當土壤含水量高于θ'D時，k值開始減小，f值降低，土壤凍融作用對土體抗剪強度的損傷增大。另外，由于土壤抗剪強度隨凍融循環次數的增加呈先減小后趨于穩定的趨勢，因此，多個水文年或多個時段計算的f值不能直接累加計算。

2 研究區與數據

2.1 研究區概況

以青藏高原東南部的高山峽谷地區為研究區（90.67°～104.2° E、28.6°～31.46° N），總面積為4.12×105km2（圖2）。研究區地勢起伏劇烈，自東向西分別為四川盆地西部、青藏高原東南緣高山峽谷地帶以及青藏高原腹地南部，海拔范圍為302～6 436 m（平均約4 000 m）。研究區中部大面積的高山峽谷區歷史構造運動強烈，地質構造發育，且巖土體破碎嚴重，是各類地質災害的高易發和高風險區，對當地生態環境安全和人類生產建設安全具有潛在威脅。

圖2 研究區位置和基本地理環境Fig. 2 Location and basic geographical environment of the study area： location （a）， terrain and distribution of observation stations （b）， soil types （c）， and land use/cover types （d）

21 世紀初，研究區多年凍土和季節凍土覆蓋全域，其中金沙江以西主要為大片島狀多年凍土，以東為季節凍土和短時凍土。受氣候變暖的影響，青藏高原總體變暖的趨勢十分明顯［51-52］。近40 年來，青藏高原升溫的速率比全球同期升溫速率高2 倍左右［2］，青藏高原東南緣是青藏高原地區升溫最顯著的區域之一［53］。因而，該地區多年凍土和季節凍土都發生了劇烈的變化和退化，對高原地表景觀形態和地表自然過程產生了深遠的影響。在此背景下，探索和查明高原寒區季節凍土的分布格局及其時空變化，是為寒區基礎自然科學研究和保障人類生產活動和生存環境安全提供科學依據的重要手段。

2.2 基礎數據

本研究收集了DEM 高程數據、土壤質地、土地利用類型等基礎資料，其中，DEM 來源于USDS EROS 數據中心（https：//srtm. csi. cgiar. org/srtmdata/）發布的90 m 空間分辨率的SRTM DEM。土壤類型空間分布數據來源于全國土壤普查辦公室編制并出版的《1∶100 萬中華人民共和國土壤圖》，基本單元為土壤亞類（研究區共有79 類，土壤深度約0.4～1.5 m），包含16 種土壤屬性數據項，如土壤深度、土壤機械組成、土壤有機質含量等，土壤屬性數據庫可在線訪問（http：//vdb3. soil. csdb. cn）。土地利用類型數據來源于中國多時相土地利用現狀數據庫。土壤類型和土地利用類型數據均下載自中國科學院資源環境科學數據中心（https：//www.resdc.cn）。FFIMS模型通過數據查找表的方式實現模型所需的土壤和植被相關參數的空間參數化。

2.3 氣象驅動數據

（1）氣溫和降水

空間分布式的氣溫和降水數據是FFIMS 模型的主要驅動參數，本研究收集了研究時段內的1 km分辨率氣溫和降水同化數據產品（a high-resolution and long-term gridded dataset for temperature and precipitation across China， HRLT）［54］，數據來源于地球和科學環境科學數據庫（https：//doi. pangaea.de）。HRLT 數據產品使用綜合統計分析方法，結合高程及其地形參數、地形濕度指數等協變量，對中國國家氣象數據共享網發布的0.5°網格逐日氣象數據參數進行空間插值，得到全國長時段高分辨率的最高氣溫、最低氣溫和降水數據產品。最高溫度的平均絕對誤差（MAE）、確定系數（R2）和Nash建模效率（NSE）分別為1.07 °C、0.98 和0.98；最低溫度的MAE、R2和NSE 分別為1.08 °C、0.99 和0.99；降水數據的MAE、R2和NSE 分別為1.30 mm、0.71和0.70。

（2）其他氣象參數

其他氣象參數包括相對濕度、潛在蒸散發、大氣壓強、風速、日照時長、0 cm 地表溫度等，數據主要來源于中國國家氣象數據共享網提供的研究區范圍內的25 個氣象站（圖2）的逐日觀測資料（http：//data.cma.cn）。采用自然鄰域空間插值方法（Natural Neighbor， NN）［55］，將氣象站點數據轉換為1 km 分辨率的空間格網數據。由于上述氣象參數具有空間大尺度的特點，利用有限的站點數據插值的方法是可行的。

2.4 站點觀測數據

對于模型精度驗證，本研究獲取了位于青藏高原中東部那曲地區附近的土壤溫/濕度監測數據集（SMTMN），該監測網絡利用ECH2O EC-TM 電磁傳感器獲得了研究區56個站點（位于本研究區的站點為16個，圖2）上的地表土壤濕度和溫度的長時間逐日地面觀測數據集，數據采集時段為2010 年8 月1 日至2016 年12 月31 日。該數據集通過利用土壤質地和土壤有機碳含量對土壤含水量和土壤溫度觀測記錄進行了校準和質量控制處理，具有較高的精度［56］。

2.5 遙感數據產品

本研究收集了兩種地表土壤水分遙感數據產品，與FFIMS 模擬的土壤含水量進行對比分析。一個是風云三號B星攜帶的微波濕度計獲取的逐日地表土壤濕度數據集（以下簡稱“FY-3B”），下載自風云衛星數據中心（http：//satellite. nsmc. org. cn），該數據集的空間分辨率為0.25°，時間范圍為2011年1月1 日至2018 年12 月31 日。另一個是融合了AMSR-E（NASA National Snow and Ice Data Center Distributed Active Archive Center）、AMSR2（Global Change Observation Mission）和SMOS（Soil Moisture and Ocean Salinity）土壤濕度數據產品的逐月地表土壤濕度數據集（以下簡稱“AAS”），下載自國家青藏高原科學數據中心（https：//data. tpdc. ac. cn），該數據集的空間分辨率為0.05°，時間范圍為2011年1月至2018年12月。

2.6 土壤溫濕度初始條件

由于模型的初始土壤水熱條件是多層的結構，而土壤溫濕度數據通常為地表單層土壤的狀態。因此，在模型輸入數據預處理模塊，預設了兩種土壤剖面一維溫度場和濕度場的計算方法。通過地表單層土壤的溫濕度數據，即可近似估算土壤垂直剖面上的溫度和水分分布。土壤剖面一維溫度場和濕度場的計算方法分別見式（9）和式（10）［5，57］。

式中：Tz為土壤深度為z時的土壤溫度（℃）；z0為地表土壤層的厚度（m），z0

3 結果與討論

3.1 凍土水熱過程數值模擬結果與驗證

FFIMS 模型模擬的部分結果展示見圖3，空間分布圖包括土壤剖面平均溫度、土壤剖面平均含水/冰量以及地表熱通量，曲線圖為逐日的區域平均土壤溫度、含水/冰量。

圖3 FFIMS模型部分模擬結果的時空分布Fig. 3 Spatiotemporal distribution of partial simulation results of the FFIMS model

陸面過程模型模擬的參量較多，而用于模型精度驗證的實測數據類型相對較少。在模型模擬過程中，各參數之間相互聯系又相互影響。一般而言，如果主要模擬參數的驗證指標良好，則通常認為模型精度是可靠的。土壤剖面溫度和含水量是凍土的關鍵特征變量，本研究主要利用青藏高原那曲地區16個野外站點的土壤溫濕度觀測數據，采用納什系數（NSE）、均方根誤差（RMSE）和決定系數（R2）指標，驗證FFIMS 模型的精度，總體精度驗證結果見圖4，逐年精度驗證結果見表2。FFIMS 模型模擬的土壤剖面溫度的納什系數和均方根誤差分別為0.89（0.87～0.91）、2.38（2.33～2.74），土壤剖面含水量的納什系數和均方根誤差分別為0.59（0.44～0.74）、0.05（0.04～0.06），該結果驗證了FFIMS 模型在刻畫凍土水熱過程中具有良好的精度。從R2指標的驗證結果看，模型模擬的地表土壤溫度和土壤含水量的時間變化趨勢在驗證時段內是十分一致的，證明了FFIMS 模型對凍土水熱過程模擬的穩定性。

表2 FFIMS模型逐年精度驗證結果Table 2 Annual precision verification results of the FFIMS model

圖4 FFIMS模型在16個凍土觀測站的精度驗證結果：土壤溫度模擬結果與實測值時間序列對比（a），土壤溫度模擬精度（b），土壤溫度模擬值與實測值的箱線圖（c），土壤含水量模擬結果與實測值時間序列對比（d），土壤含水量模擬精度（e），土壤含水量模擬值與實測值的箱線圖（f）Fig. 4 Overall accuracy verification results of the FFIMS model： comparison of time series between simulated soil temperature results and measured values （a）， simulation accuracy of soil temperature （b）， boxplot of simulated and observed soil temperatures （c）， comparison of time series between simulated soil water content results and measured values （d），simulation accuracy of soil water content （e）， and boxplot of simulated and observed soil water contents （f）

圖5展示了FFIMS模型模擬的土壤表層水含量與現有數據產品在逐日和逐月尺度上的對比驗證結果?？梢园l現，模型模擬土壤水含量的變化趨勢與現有的土壤濕度遙感數據產品具有較好的一致性。但是，與逐日的遙感數據產品的劇烈變化幅度相比，模型模擬的地表土壤水含量表現出了更穩定的變化趨勢。

圖5 FFIMS模型模擬的土壤表層水含量與遙感數據產品對比驗證結果Fig. 5 Comparison of surface soil water content simulated by the FFIMS model and remote sensing data products：daily FY-3B soil moisture data （a）， and monthly AAS soil moisture data （b）

3.2 凍土變化特征分析

（1）時間變化特征

本研究主要通過土壤剖面平均溫度、土壤剖面平均含水/冰量和地表熱通量等參數描述凍土系統水熱過程，圖6 展示了以上多種參量的逐日和逐年變化曲線，圖中灰色區域的上下邊界分別表示該參量在研究區范圍內的最大和最小值，紅色折線表示年均值。其中，地表熱通量指地表土壤單位面積的熱量交換（W·m-2），具有方向性，正值表示大氣向土壤傳遞熱量。

圖6 凍土水熱參量逐日和逐年時間變化趨勢Fig. 6 Trend of daily and annual changes in the water and heat parameters of frozen soil： profile temperature （a），liquid water content （b）， ice content （c）， and surface heat flux （d）

在本研究時段內，研究區氣溫增幅明顯，年平均氣溫自4.0 ℃上升到5.9 ℃，平均增幅為0.2 ℃·a-1?？偟膩碚f，從逐日的凍融過程看，地表土壤凍融循環中的各水熱參量呈周期性的規律變化，逐年變化幅度不明顯。從研究區平均的角度看，近十年，青藏高原東南緣凍土系統呈現明顯的消退趨勢，主要表現為土壤剖面溫度上升、土壤剖面含水量增加等。其中，土壤剖面的年均溫度從3.5 ℃上升到6.2 ℃，平均溫度上升率達1.9 ℃·（10a）-1；土壤剖面含水量從0.2 m3·m-3上升到0.26 m3·m-3，平均上升率為0.07 m3·m-3·（10a）-1；土壤剖面含冰量呈現波動降低趨勢，但變化幅度較低，總體趨勢平穩；地表熱通量則呈顯著上升趨勢，從2.9 W·m-2上升到12.0 W·m-2，平均上升率為6.8 W·m-2·（10a）-1。地表土壤熱通量上升，說明大氣向陸地輸送的熱量增加，且熱通量數值為正，意味著土壤溫度總體上升的趨勢仍將持續。

（2）趨勢率空間分布特征

空間傾向率方法是在普通線性傾向率方法的基礎上，擴展至空間網格尺度上的方法，即在每一個空間柵格上建立相關變量的時間序列，對時間序列逐一進行線性擬合，獲取線性擬合的斜率（趨勢率）與顯著性檢驗值，從而形成空間分布的趨勢特征。本研究利用空間傾向率方法分析研究區凍土相關特征參數的時空演化規律，重點分析各水熱參量變化的空間分布特征。根據逐年土壤剖面平均溫度、含水/冰量、地表熱通量空間分布數據進行空間傾向率計算，獲得其時間變化趨勢的空間分布，結果如圖7所示。

從凍土水熱參量變化趨勢率的空間分布上看，各種參量的空間分布具有顯著的異質性，空間分布特征受地形地貌、土壤質地等地表景觀因素以及氣候要素的影響。由于土壤溫度是其他參量變化的主導因素，土壤溫度與土壤含水/冰量等要素的趨勢率在空間上的分布特征具有一定程度的相似性。其中，土壤溫度變化趨勢的空間分區特征最為明顯，研究區西部地處青藏高原腹地的區域，除南部低洼地區，土壤溫度總體呈顯著上升的趨勢；中部川藏高山峽谷地區地形多變，土壤溫度變化的空間特征復雜；東部四川盆地的土壤溫度呈不顯著下降趨勢或基本無變化，但與四川盆地交界的青藏高原東南邊緣地帶，是研究區內土壤溫度上升幅度最大的區域。土壤含水量總體上呈現上升趨勢，研究區大部分區域為顯著上升趨勢，局部為不顯著上升，研究區南部洼地小面積的區域存在下降趨勢。土壤含冰量總體呈現不顯著下降的變化趨勢，在空間分布上與土壤含水量上升趨勢的分布特征較一致，其中呈顯著下降趨勢的地區多位于研究區東部的地勢較低的溝谷地帶。隨著年均氣溫的升高，研究區地表土壤熱通量呈逐年增加的趨勢，且大部分地區較為顯著，局部地區呈現下降的趨勢。由于地表熱通量隨時間變化較劇烈，即使在同一天也可能出現較大的波動，因此其變化趨勢的空間分布無明顯規律性特征。

3.3 土壤凍融作用表征指數

根據FFIMS 模型模擬的精細凍土水熱過程參量，按照本研究提出的凍融作用參數化表征方法，計算了研究區土體抗剪強度損傷系數。圖8展示了不同凍融循環次數（1 個水文年為1 次凍融循環，圖中以n表示）下的土體抗剪強度損傷系數空間分布。

圖8 不同凍融循環次數下的土體抗剪強度損傷系數空間分布Fig. 8 Spatial distribution of the damage coefficient of soil shear strength under different freeze-thaw cycles

基于凍土凍融循環水熱參量計算的土體抗剪強度的損傷系數，主要反映了區域地表土壤凍融作用對土體抗剪強度的相對損傷程度。損傷系數越低，表示損傷程度越大，邊坡失穩發生凍融災害的風險或概率（相對）加大。從時間變化上看，隨著凍融循環次數的增加，凍融作用的損傷程度不斷加劇，但其空間分布特征總體較穩定。從空間分布上看，土體抗剪強度損傷系數具有較強的地帶性分布規律。高值區域集中分布在研究區東部和南部低洼峽谷地區，這些地區凍結天數較少甚至未發生凍結，土壤凍融作用相對較弱。無論在研究區西部的青藏高原腹地、中部的川藏高山峽谷地區，還是研究區東部低洼地區，均存在損傷程度較高（低值）的區域。其中，青藏高原東部邊緣與四川盆地接壤地帶有成片低值區域分布。受青藏高原東部邊緣多年凍土退化程度較大的影響，該區域地表土壤凍融作用整體較強。

本研究以凍融循環5 次的計算結果為例，進一步結合地形、植被類型進行了統計分析（圖9）。統計結果表明，地表土壤凍融作用對土體的損傷程度與高程有較顯著的相關性，高程低于1 000 m 的區域損傷程度基本較小，而青藏高原腹地和高山峽谷區高程在5 000 m 左右的區域，土體強度受凍融作用影響最大。由于研究區內的人類活動多處于地勢低洼的地區，因此土地利用類型為農田和建設用地的區域，地表土壤凍融作用的損傷程度較小。反之，高海拔山區的草地和林地，尤其是地表裸露的區域（鹽堿地、戈壁、裸土/石等），土體抗剪強度在凍融循環的作用下損傷程度較大。另外，通過對比不同植被覆蓋率下的損傷系數發現，植被覆蓋度越高，損傷系數越高，土壤凍融作用的損傷程度越低。

圖10 展示了研究區平均土體抗剪強度和土壤冰融化量隨凍融循環次數增加的變化。其中，年均損傷系數和累積損傷系數的區別在于，計算累積損傷系數的基準年份是研究時段的起始年份。通過定量統計可以發現，累積損傷系數隨凍融循環次數的增加而降低，但隨凍融循環次數的增加，累積損傷系數降低的程度逐漸減小，直至平穩的趨勢。這表明凍融循環導致土體損傷失穩的作用具有一定的限度，不會成為引起凍融災害發生的主導因素。

圖10 土體抗剪強度損傷系數和土壤冰融量隨時間的變化Fig. 10 Variations of the damage coefficient of soil shear strength and melted soil ice with time

另外，需要特別注意的是，本研究計算得到的損傷系數并不能反映土體抗剪強度本身大小或邊坡固有穩定性特點。例如，損傷系數較低的區域，表示凍融作用對土體的損傷程度大，但并非意味著邊坡穩定性差或凍融災害發生概率高，僅表示該區域地表巖土體受到的凍土凍融作用更強烈。除了凍融作用外，土體本身的抗剪強度對邊坡穩定性起到決定性作用，凍融循環起到了弱化其相對強度的作用。

4 結論

針對當前分布式凍土過程模擬研究的不足，彌補區域尺度凍土研究在自然災害風險評估與防控領域中應用研究的短板，本研究以青藏高原東南緣的高山峽谷區及其周邊地區為研究區，重點關注季節凍土凍融循環及其水熱傳輸過程的系統性表達，建立了適用于青藏高原高海拔凍土區的空間全分布式的凍土水熱耦合過程數值模型（FFIMS 模型）。在獲取高精度的凍土系統分布式水熱過程參量的基礎上，提出凍土凍融作用的空間參數化方法，得到以下主要結論：

（1）基于物理機制的FFIMS 模型對地表凍融循環中的凍土水熱過程的模擬效果較好，凍土系統主要特征參量（土壤溫度和土壤含水量）經過驗證精度良好，可以為精細刻畫空間分布式凍土過程提供有效的模型方法。

（2）隨著氣溫升高，青藏高原東南緣季節凍土變化劇烈，空間異質性較強，季節凍土除了周期性凍融循環外，總體呈退化的趨勢，表現為土壤溫度顯著上升、土壤含水量顯著上升、土壤含冰量不顯著下降等，為凍融作用下的巖土體結構的抗剪強度變化增加了更多不確定性因素。

（3）另外，還進行了地表土壤凍融作用空間參數化表征方法研究，創新性地提出了土體抗剪強度損傷系數，通過統計分析表明該指標在數值表達和空間分布等方面是合理可行的。利用高精度的凍土過程特征參量，可以有效進行凍土凍融作用的參數化表達。

本研究成果可以為凍土相關研究提供新的思路，同時，為寒區區域尺度下的相關基礎研究和災害風險評估與災害防治等領域的研究提供凍土系統動態演化的數據和技術支撐。

致謝：本文使用的土壤溫濕度觀測數據和遙感反演的土壤水分數據下載自國家青藏高原科學數據中心（http://data.tpdc.ac.cn），土壤類型和土地利用類型數據下載自中國科學院資源環境科學數據中心（https://www.resdc.cn），氣象驅動數據下載自地球和環境科學數據庫（https://doi.pangaea.de）以及國家氣象數據共享網（http://data.cma.cn），在此一并表示感謝。