高原寒區(qū)地表土壤凍融作用的空間參數(shù)化表征研究

高會(huì)然， 許 沖， 張萬(wàn)昌， 易亞寧， 肖子亢

（1. 應(yīng)急管理部 國(guó)家自然災(zāi)害防治研究院，北京 100085； 2. 中國(guó)科學(xué)院 空天信息創(chuàng)新研究院，北京 100094；3. 復(fù)合鏈生自然災(zāi)害動(dòng)力學(xué)應(yīng)急管理部重點(diǎn)實(shí)驗(yàn)室，北京 100085）

0 引言

在全球氣候持續(xù)變暖的背景下，冰凍圈領(lǐng)域已成為研究熱點(diǎn)［1-4］。當(dāng)前，冰凍圈研究不僅關(guān)注其在氣候系統(tǒng)中的作用，還更加關(guān)注對(duì)人類社會(huì)可持續(xù)發(fā)展的現(xiàn)實(shí)影響和潛在威脅［5-7］。作為冰凍圈的重要組成部分，凍土的變化或退化對(duì)生態(tài)環(huán)境系統(tǒng)、人類生活和生產(chǎn)安全等方面的直接或間接風(fēng)險(xiǎn)，這些研究方向也越來(lái)越受到學(xué)者的關(guān)注［8-9］。

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

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

當(dāng)前，基于物理機(jī)理的凍土水熱傳輸過(guò)程數(shù)值模型研究取得了較大的進(jìn)展，部分陸面過(guò)程模型通過(guò)模型集成的方式，具備了精細(xì)刻畫地表土壤凍融循環(huán)多參量和多過(guò)程的能力［39-40］。但是，目前凍土變化監(jiān)測(cè)與模擬研究成果在陸面環(huán)境要素響應(yīng)和凍土災(zāi)害防治等領(lǐng)域的應(yīng)用能力仍處于較低的水平。一方面，高時(shí)空分辨率凍土特征參量數(shù)據(jù)獲取難度相對(duì)較大，限制了大范圍長(zhǎng)時(shí)段凍土凍融循環(huán)及其凍融作用的定量分析；另一方面，地表土壤凍融作用作為凍土水熱過(guò)程與陸面景觀格局演變相互聯(lián)系的重要紐帶，其表征指標(biāo)識(shí)別與空間參數(shù)化的定量研究不足。區(qū)域尺度凍土凍融作用的參數(shù)化表征及其進(jìn)一步在災(zāi)害風(fēng)險(xiǎn)評(píng)估與防控領(lǐng)域中的應(yīng)用研究仍是空白。

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

1 理論與方法

1.1 總體方法框架

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

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

1.2 分布式凍土水熱過(guò)程數(shù)值建模

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

（1） 植被冠層

植被及其冠層設(shè)計(jì)為單層單節(jié)點(diǎn)結(jié)構(gòu)，同時(shí)考慮了植被生長(zhǎng)和季節(jié)變化特征，冠層能量與水平衡方程如式（1）和式（2）所示。

式中：ρa(bǔ)、ca和T分別為大氣密度（kg·m-3）、大氣比熱容（J·kg-1·℃-1）和冠層溫度（℃）；ke為冠層熱量傳導(dǎo)系數(shù)（m2·s-1）；Hl為冠層向大氣傳輸?shù)臒崃浚╓·m-3）；ρv為蒸氣密度（kg·m-3）；El為冠層蒸散發(fā)量（mm）。

（2） 積雪

將積雪設(shè)計(jì)為單層雙節(jié)點(diǎn)的結(jié)構(gòu)，考慮積雪的累積和升華以及水熱在雪層中的傳導(dǎo)過(guò)程，積雪能量平衡方程可用式（3）表示。

式中：ρsp、wsp和ksp分別為雪密度（kg·m-3）、體積含水量（m3·m-3）和熱傳導(dǎo)系數(shù)（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） 土壤

設(shè)計(jì)為不同網(wǎng)格不同深度的多層土壤結(jié)構(gòu)，考慮了凍融過(guò)程對(duì)水熱通量的影響，土壤能量和水平衡方程可分別表示為式（4）和式（5）。

式中：Cs為土壤比熱容（J·kg-1·℃-1）；ks為土壤熱傳導(dǎo)系數(shù)（W·m-1·℃-1）；K為土壤導(dǎo)水率（m·s-1）；ψ為土壤水勢(shì)（m）；U為土壤水通量的源/匯項(xiàng)（m3·m-3·s-1）。

在FFIMS 模型的水熱平衡方程的求解過(guò)程中，將植被冠層、積雪和土壤剖面分為有限層數(shù)，每層用一個(gè)節(jié)點(diǎn)表示，每個(gè)節(jié)點(diǎn)中的各水熱分量的存儲(chǔ)量以各節(jié)點(diǎn)所在層的土壤厚度為準(zhǔn)。利用隱式有限差分方程組求解每個(gè)節(jié)點(diǎn)以及節(jié)點(diǎn)間的能量和水的平衡方程［43］。模型中各個(gè)子過(guò)程的理論和方法涉及大量參數(shù)和變量，包括一些可調(diào)參數(shù)、經(jīng)驗(yàn)常數(shù)和物理常數(shù)等。其中，可調(diào)參數(shù)主要是指受研究區(qū)域環(huán)境條件影響的參數(shù)，需要根據(jù)研究區(qū)地形、氣象、生態(tài)等條件調(diào)整。本研究建立的模型中主要的參數(shù)取值及其說(shuō)明見表1。

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

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

1.3 凍融作用參數(shù)化表征方法

地表土壤的凍融作用直接影響或引起土體分裂以及土壤顆粒的重新組合，從而影響土體的抗剪強(qiáng)度，導(dǎo)致土體從穩(wěn)定狀態(tài)轉(zhuǎn)向相對(duì)不穩(wěn)定狀態(tài)，威脅地區(qū)的邊坡穩(wěn)定、工程安全等［44-45］。但是，凍融作用對(duì)土抗剪強(qiáng)度以及相關(guān)指標(biāo)的影響程度目前仍無(wú)定論。根據(jù)本研究獲取的精細(xì)凍土特征參量，結(jié)合現(xiàn)有的理論研究基礎(chǔ)，可以實(shí)現(xiàn)定量化表征凍土凍融作用對(duì)土體抗剪強(qiáng)度指標(biāo)的相對(duì)損傷程度。因此，本研究以土壤剖面含水量和土壤剖面含冰量?jī)蓚€(gè)關(guān)鍵凍土特征參量為主要參數(shù)，構(gòu)建土體抗剪強(qiáng)度的損傷系數(shù)。

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

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

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

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

2 研究區(qū)與數(shù)據(jù)

2.1 研究區(qū)概況

以青藏高原東南部的高山峽谷地區(qū)為研究區(qū)（90.67°～104.2° E、28.6°～31.46° N），總面積為4.12×105km2（圖2）。研究區(qū)地勢(shì)起伏劇烈，自東向西分別為四川盆地西部、青藏高原東南緣高山峽谷地帶以及青藏高原腹地南部，海拔范圍為302～6 436 m（平均約4 000 m）。研究區(qū)中部大面積的高山峽谷區(qū)歷史構(gòu)造運(yùn)動(dòng)強(qiáng)烈，地質(zhì)構(gòu)造發(fā)育，且?guī)r土體破碎嚴(yán)重，是各類地質(zhì)災(zāi)害的高易發(fā)和高風(fēng)險(xiǎn)區(qū)，對(duì)當(dāng)?shù)厣鷳B(tài)環(huán)境安全和人類生產(chǎn)建設(shè)安全具有潛在威脅。

圖2 研究區(qū)位置和基本地理環(huán)境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 世紀(jì)初，研究區(qū)多年凍土和季節(jié)凍土覆蓋全域，其中金沙江以西主要為大片島狀多年凍土，以東為季節(jié)凍土和短時(shí)凍土。受氣候變暖的影響，青藏高原總體變暖的趨勢(shì)十分明顯［51-52］。近40 年來(lái)，青藏高原升溫的速率比全球同期升溫速率高2 倍左右［2］，青藏高原東南緣是青藏高原地區(qū)升溫最顯著的區(qū)域之一［53］。因而，該地區(qū)多年凍土和季節(jié)凍土都發(fā)生了劇烈的變化和退化，對(duì)高原地表景觀形態(tài)和地表自然過(guò)程產(chǎn)生了深遠(yuǎn)的影響。在此背景下，探索和查明高原寒區(qū)季節(jié)凍土的分布格局及其時(shí)空變化，是為寒區(qū)基礎(chǔ)自然科學(xué)研究和保障人類生產(chǎn)活動(dòng)和生存環(huán)境安全提供科學(xué)依據(jù)的重要手段。

2.2 基礎(chǔ)數(shù)據(jù)

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

2.3 氣象驅(qū)動(dòng)數(shù)據(jù)

（1）氣溫和降水

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

（2）其他氣象參數(shù)

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

2.4 站點(diǎn)觀測(cè)數(shù)據(jù)

對(duì)于模型精度驗(yàn)證，本研究獲取了位于青藏高原中東部那曲地區(qū)附近的土壤溫/濕度監(jiān)測(cè)數(shù)據(jù)集（SMTMN），該監(jiān)測(cè)網(wǎng)絡(luò)利用ECH2O EC-TM 電磁傳感器獲得了研究區(qū)56個(gè)站點(diǎn)（位于本研究區(qū)的站點(diǎn)為16個(gè)，圖2）上的地表土壤濕度和溫度的長(zhǎng)時(shí)間逐日地面觀測(cè)數(shù)據(jù)集，數(shù)據(jù)采集時(shí)段為2010 年8 月1 日至2016 年12 月31 日。該數(shù)據(jù)集通過(guò)利用土壤質(zhì)地和土壤有機(jī)碳含量對(duì)土壤含水量和土壤溫度觀測(cè)記錄進(jìn)行了校準(zhǔn)和質(zhì)量控制處理，具有較高的精度［56］。

2.5 遙感數(shù)據(jù)產(chǎn)品

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

2.6 土壤溫濕度初始條件

由于模型的初始土壤水熱條件是多層的結(jié)構(gòu)，而土壤溫濕度數(shù)據(jù)通常為地表單層土壤的狀態(tài)。因此，在模型輸入數(shù)據(jù)預(yù)處理模塊，預(yù)設(shè)了兩種土壤剖面一維溫度場(chǎng)和濕度場(chǎng)的計(jì)算方法。通過(guò)地表單層土壤的溫濕度數(shù)據(jù)，即可近似估算土壤垂直剖面上的溫度和水分分布。土壤剖面一維溫度場(chǎng)和濕度場(chǎng)的計(jì)算方法分別見式（9）和式（10）［5，57］。

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

3 結(jié)果與討論

3.1 凍土水熱過(guò)程數(shù)值模擬結(jié)果與驗(yàn)證

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

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

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

表2 FFIMS模型逐年精度驗(yàn)證結(jié)果Table 2 Annual precision verification results of the FFIMS model

圖4 FFIMS模型在16個(gè)凍土觀測(cè)站的精度驗(yàn)證結(jié)果：土壤溫度模擬結(jié)果與實(shí)測(cè)值時(shí)間序列對(duì)比（a），土壤溫度模擬精度（b），土壤溫度模擬值與實(shí)測(cè)值的箱線圖（c），土壤含水量模擬結(jié)果與實(shí)測(cè)值時(shí)間序列對(duì)比（d），土壤含水量模擬精度（e），土壤含水量模擬值與實(shí)測(cè)值的箱線圖（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模型模擬的土壤表層水含量與現(xiàn)有數(shù)據(jù)產(chǎn)品在逐日和逐月尺度上的對(duì)比驗(yàn)證結(jié)果?？梢园l(fā)現(xiàn)，模型模擬土壤水含量的變化趨勢(shì)與現(xiàn)有的土壤濕度遙感數(shù)據(jù)產(chǎn)品具有較好的一致性。但是，與逐日的遙感數(shù)據(jù)產(chǎn)品的劇烈變化幅度相比，模型模擬的地表土壤水含量表現(xiàn)出了更穩(wěn)定的變化趨勢(shì)。

圖5 FFIMS模型模擬的土壤表層水含量與遙感數(shù)據(jù)產(chǎn)品對(duì)比驗(yàn)證結(jié)果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）時(shí)間變化特征

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

圖6 凍土水熱參量逐日和逐年時(shí)間變化趨勢(shì)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）

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

（2）趨勢(shì)率空間分布特征

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

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

3.3 土壤凍融作用表征指數(shù)

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

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

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

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

圖10 展示了研究區(qū)平均土體抗剪強(qiáng)度和土壤冰融化量隨凍融循環(huán)次數(shù)增加的變化。其中，年均損傷系數(shù)和累積損傷系數(shù)的區(qū)別在于，計(jì)算累積損傷系數(shù)的基準(zhǔn)年份是研究時(shí)段的起始年份。通過(guò)定量統(tǒng)計(jì)可以發(fā)現(xiàn)，累積損傷系數(shù)隨凍融循環(huán)次數(shù)的增加而降低，但隨凍融循環(huán)次數(shù)的增加，累積損傷系數(shù)降低的程度逐漸減小，直至平穩(wěn)的趨勢(shì)。這表明凍融循環(huán)導(dǎo)致土體損傷失穩(wěn)的作用具有一定的限度，不會(huì)成為引起凍融災(zāi)害發(fā)生的主導(dǎo)因素。

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

另外，需要特別注意的是，本研究計(jì)算得到的損傷系數(shù)并不能反映土體抗剪強(qiáng)度本身大小或邊坡固有穩(wěn)定性特點(diǎn)。例如，損傷系數(shù)較低的區(qū)域，表示凍融作用對(duì)土體的損傷程度大，但并非意味著邊坡穩(wěn)定性差或凍融災(zāi)害發(fā)生概率高，僅表示該區(qū)域地表巖土體受到的凍土凍融作用更強(qiáng)烈。除了凍融作用外，土體本身的抗剪強(qiáng)度對(duì)邊坡穩(wěn)定性起到?jīng)Q定性作用，凍融循環(huán)起到了弱化其相對(duì)強(qiáng)度的作用。

4 結(jié)論

針對(duì)當(dāng)前分布式凍土過(guò)程模擬研究的不足，彌補(bǔ)區(qū)域尺度凍土研究在自然災(zāi)害風(fēng)險(xiǎn)評(píng)估與防控領(lǐng)域中應(yīng)用研究的短板，本研究以青藏高原東南緣的高山峽谷區(qū)及其周邊地區(qū)為研究區(qū)，重點(diǎn)關(guān)注季節(jié)凍土凍融循環(huán)及其水熱傳輸過(guò)程的系統(tǒng)性表達(dá)，建立了適用于青藏高原高海拔凍土區(qū)的空間全分布式的凍土水熱耦合過(guò)程數(shù)值模型（FFIMS 模型）。在獲取高精度的凍土系統(tǒng)分布式水熱過(guò)程參量的基礎(chǔ)上，提出凍土凍融作用的空間參數(shù)化方法，得到以下主要結(jié)論：

（1）基于物理機(jī)制的FFIMS 模型對(duì)地表凍融循環(huán)中的凍土水熱過(guò)程的模擬效果較好，凍土系統(tǒng)主要特征參量（土壤溫度和土壤含水量）經(jīng)過(guò)驗(yàn)證精度良好，可以為精細(xì)刻畫空間分布式凍土過(guò)程提供有效的模型方法。

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

（3）另外，還進(jìn)行了地表土壤凍融作用空間參數(shù)化表征方法研究，創(chuàng)新性地提出了土體抗剪強(qiáng)度損傷系數(shù)，通過(guò)統(tǒng)計(jì)分析表明該指標(biāo)在數(shù)值表達(dá)和空間分布等方面是合理可行的。利用高精度的凍土過(guò)程特征參量，可以有效進(jìn)行凍土凍融作用的參數(shù)化表達(dá)。

本研究成果可以為凍土相關(guān)研究提供新的思路，同時(shí)，為寒區(qū)區(qū)域尺度下的相關(guān)基礎(chǔ)研究和災(zāi)害風(fēng)險(xiǎn)評(píng)估與災(zāi)害防治等領(lǐng)域的研究提供凍土系統(tǒng)動(dòng)態(tài)演化的數(shù)據(jù)和技術(shù)支撐。

致謝：本文使用的土壤溫濕度觀測(cè)數(shù)據(jù)和遙感反演的土壤水分?jǐn)?shù)據(jù)下載自國(guó)家青藏高原科學(xué)數(shù)據(jù)中心（http://data.tpdc.ac.cn），土壤類型和土地利用類型數(shù)據(jù)下載自中國(guó)科學(xué)院資源環(huán)境科學(xué)數(shù)據(jù)中心（https://www.resdc.cn），氣象驅(qū)動(dòng)數(shù)據(jù)下載自地球和環(huán)境科學(xué)數(shù)據(jù)庫(kù)（https://doi.pangaea.de）以及國(guó)家氣象數(shù)據(jù)共享網(wǎng)（http://data.cma.cn），在此一并表示感謝。