李穎，吳靜，李純斌，秦格霞

（甘肅農(nóng)業(yè)大學資源與環(huán)境學院，甘肅 蘭州 730070）

青藏高原（以下簡稱高原）占中國陸地面積的1/4，海拔均在4000 m以上，是世界上海拔最高面積最大的高原［1］。在全球氣候變暖的大背景下，高原的地表能量研究成為熱點課題［2-4］。地表層土壤熱通量（ground surface soil heat flux，G0）是地表能量平衡（surface energy balance，SEB）的重要參數(shù)之一，它可以解釋地表淺層與深層之間傳輸?shù)哪芰浚?］。同時對地氣間的物質(zhì)能量傳輸及分配有重要影響［6］。此外，G0也是地表蒸散發(fā)模擬的重要因子，會極大地影響地表蒸散發(fā)的模擬精度［7］。一直以來對高原G0的研究從未中斷，高原腹地的G0主要受到地氣溫差的影響，同時下墊面狀況也起到重要作用［8］；利用超高分辨率掃描輻射計（advanced very high resolution radiometer，AVHRR）和中分辨率成像光譜儀（moderate resolution image spectroradiometer，MODIS）反演地表特征參數(shù)和地表通量，發(fā)現(xiàn)G0空間分布與地表特征參數(shù)（植被指數(shù)、植被覆蓋度和地表反照率）的空間分布有很好的對應關系［9］。由此可見，高原G0與下墊面有非常密切的關系。高原分布著占全國草地總面積38%的天然草地，不僅具有防風固沙、涵養(yǎng)水源、調(diào)節(jié)碳循環(huán)及氣候變化、維護生物多樣性等多種生態(tài)服務功能［10］，同時還是當?shù)鼐用裆姘l(fā)展和畜牧業(yè)發(fā)展的基礎［11］。通過研究高原草地G0，不僅有助于加深對高原G0的理解，而且對高原土壤-植被-大氣系統(tǒng)的能量平衡研究有所幫助。

目前對高原地表層土壤熱通量G0的研究主要有兩種手段，一是基于站點觀測數(shù)據(jù)估算分析G0，二是結合遙感影像數(shù)據(jù)分析時空尺度較大的G0，中心思想是通過獲取G0與凈輻射（net radiation，Rn）的比率，以此估計G0。利用不同的地表參數(shù)估算G0/Rn，如基于植被參數(shù)的模型，有歸一化植被指數(shù)（normalized difference vegetation index，NDVI）和葉面積指數(shù)（leaf area index，LAI）模型［12］；除此之外考慮地表溫度和地表反照率的模型，如地表能量平衡模型（surface energy balance model，SEBAL）［13］和Ma模型［14］。有學者對這些模型做了對比分析，對于特殊的高原氣候條件，反演地表層土壤熱通量最適用的方法是Ma模型［15］。由于高原的地理氣候條件，實測數(shù)據(jù)難以獲取，因此利用遙感數(shù)據(jù)估算G0是非常有效的手段。本研究利用Ma模型反演高原2003-2018年地表層土壤熱通量（G0），分析2003-2018年高原G0的時空變化特征，探討不同草地類型的G0年際變化以及季節(jié)變化。

1 材料與方法

1.1 研究區(qū)概況

青藏高原（26°00′12″-39°46′50″N，73°18′52″-104°46′59″E）總體地勢西高東低，地形復雜，地勢險峻多變。草地是高原主要的生態(tài)系統(tǒng)類型，分布情況如圖1所示，主要以高寒草原類和高寒草甸類為主，占高原草地面積的73.8%。

1.2 數(shù)據(jù)來源及處理

1.2.1 草地類型數(shù)據(jù) 草地類型數(shù)據(jù)來自蘇大學［16］團隊編制的1∶1000000的中國草地類型數(shù)據(jù)，利用高原行政區(qū)域提取高原的草地類型，分類合并，得到高原的草地類型圖（圖1）。

圖1 青藏高原草地類型及觀測點分布Fig.1 Grassland types and distribution of observation points on the Qinghai Tibet Plateau

1.2.2 站點數(shù)據(jù) 選擇7個高原觀測土壤熱通量的站點（圖1和表1），其中大沙龍站，阿柔站，峨堡站，景陽嶺站，埡口站數(shù)據(jù)來源于國家青藏高原科學數(shù)據(jù)中心（http：//data.tpdc.ac.cn/zh-hans/），本數(shù)據(jù)由“黑河生態(tài)水文遙感試驗（HiWATER）”產(chǎn)生［17-18］；西大灘站，唐古拉站數(shù)據(jù)來源于國家冰川凍土沙漠科學數(shù)據(jù)中心（http：//www.crensed.ac.cn/portal/）［19-20］。唐古拉站和西大灘站的地表溫度（temperature，T）用史蒂芬-玻爾茲曼公式計算獲得［21］。

表1 觀測站點信息Table 1 The information of observation sites

式中：T為地表溫度（單位：K）；Rul和Rdl分別為地表向上長波輻射和地表向下長波輻射（單位：W·m-2）；εg為地表比輻射率，取值為0.96；σ為Stefan-Boltzmann常數(shù)，取5.67×10-8W·m-2·K-4。

1.2.3 遙感數(shù)據(jù) 地表層土壤熱通量的遙感反演數(shù)據(jù)包括中國西部1 km全天候地表溫度數(shù)據(jù)V1，覆蓋整個青藏高原，時間分辨率為逐日，空間分辨率為1 km［22］；MODIS數(shù)據(jù)來自于美國航空航天局，包括MOD09GA和MOD13Q1。其中MOD09GA提供1～7波段的表面光譜反射率估計值，主要用于地表反照率的計算；MOD13Q1提供NDVI、紅光和近紅外波段地表反射率值，NDVI用于地表發(fā)射率的反演，紅光和近紅外波段的地表反射率用于反演區(qū)域改良的土壤調(diào)整植被指數(shù)（modified soil adjusted vegetation index，MSAVI）；中國區(qū)域地面氣象要素驅動數(shù)據(jù)集（1979-2018）來自國家青藏高原科學數(shù)據(jù)中心，該數(shù)據(jù)集主要提供近地面氣溫、近地面氣壓、近地面空氣比濕、近地面全風速、地面向下短波輻射、地面向下長波輻射和地面降水率共7個要素［23-26］。所需遙感數(shù)據(jù)具體介紹如表2所示。

表2 遙感數(shù)據(jù)信息Table 2 The information of remote sensing data

1.3 研究方法

首先利用單點實測數(shù)據(jù)，計算單點地表土壤熱通量；其次根據(jù)遙感反演模型，利用遙感影像數(shù)據(jù)反演得到高原地表層土壤熱通量，分析不同草地類型的地表層土壤熱通量，具體流程如圖2所示。

圖2 技術路線Fig.2 Technology road

1.3.1 單點實測數(shù)據(jù)計算方法 采用以下公式計算站點地表層土壤熱通量：

式中：t為時間（單位：s）；z為土壤深度（向下為正，單位：m）；T為土壤溫度（單位：K）；Cs為土壤比熱容（單位：J·kg-1·K-1）；ρs為土壤密度（單位：kg·m-3）；λs為土壤熱傳導系數(shù)（單位：W·K-1·m-1）；G為土壤熱通量（向下為正，單位：W·m-2）［27］。

結合公式2）與3），且考慮青藏高原近地表存在土壤凍融現(xiàn)象，因此利用5 cm的土壤熱通量計算地表層土壤熱通量（G0）［28］：

式中：zref為參考深度（基于本研究的實測數(shù)據(jù)，不同站點選擇不同的參考深度）；G(zref)為參考深度處的土壤熱通量；ρw為液態(tài)水的密度，為1.0×103kg·m-3；Lf為水的凍結-融化潛熱常數(shù)，為3.337×105J·kg-1；θw為體積含水量；Tav為5 cm土壤溫度與地表溫度的平均值。在土壤的整個凍融過程中，對不同凍融階段進行不同的計算［29］：

1）當5 cm土壤層處于完全融化階段時，土壤體積熱容量計算公式如下：

式中：ρdryCdry為干土的熱容量，約為0.90×106J·m-3·K-1［30］；ρwCw為液態(tài)水的熱容量，約為4.2×106J·m-3·K-1；θw5cm為深度5 cm處未凍水體積含水量（單位：m3·m-3）。

2）當5 cm深度以上的土壤層處于凍融循環(huán)階段時，此時土壤體積熱容量為：

式中：ρiCi為冰的熱容量，約為1.89×106J·m-3·K-1；θi5cm為5 cm處的土壤含冰量。在凍融循環(huán)階段，土壤中的含冰量可以根據(jù)對應深度的土壤中相鄰時刻之間未凍水含量來計算［31］，公式如下：

3）當土層5 cm以上完全凍結時，可以將土壤含水量和含冰量視為常數(shù)［32］。土壤中的含冰量由凍結開始前和結束后含水量的差值求得，其土壤體積熱容量仍用公式6）進行計算。

1.3.2 遙感反演方法 Ma總模型公式如下：

式中：Ts為地表溫度（單位：℃）；α是地表反照率；MSAVI為改良的土壤調(diào)整植被指數(shù)［15］；Rn是凈輻射（單位：W·m-2）。

地表反照率α基于MOD09GA的波段反射率數(shù)據(jù)計算得出［33］：

式中：R1、R2、R3、R4、R5和R7為MOD09GA的第1、2、3、4、5和7波段的地表反射率。

MSAVI通過以下公式計算得出［34］：

式中：RED和NIR分別為MOD13Q1的紅光波段和近紅外波段的地表反射率值。

凈輻射Rn可以通過以下公式計算獲?。?5］：

式中：向下短波輻射（downward short wave radiation，DSR）和向下長波輻射（downward long wave radiation，DLR）分別是向下短波輻射和向下長波輻射（單位：W·m-2）；εs是地表發(fā)射率，根據(jù)以下公式計算［36］：

式中：NDVI為歸一化植被指數(shù)，通過植被指數(shù)產(chǎn)品MOD13Q1提?。籪c為植被蓋度，可以根據(jù)NDVI估算，公式如下［37］：

式中：NDVImax和NDVImin分別是完全植被覆蓋和裸露土壤的NDVI值。

2 結果與分析

2.1 站點地表層土壤熱通量

由于G0無法直接觀測得到，因此本研究基于站點淺層的土壤熱通量觀測值G檢驗站點G0的計算結果。圖3a是唐古拉站點2015年G0與5 cm處的土壤熱通量G5之間的散點分布關系，各站點G0與淺層土壤熱通量G的線性關系式以及R2的值如表3所示，除阿柔站，其余站點的G0與G無截距線性關系，斜率略大于1，說明G0略大于G5，兩者之間的R2均大于0.8。圖3b是唐古拉站點G0與G52015年年均日變化曲線，兩者整體變化曲線呈現(xiàn)倒“U”形狀，在夜晚相較于白天變化較為平緩。在7：00之前，G5略大于G0，并且G0比G5首先達到最大值。以上結論與之前的學者研究結果一致［38］，可以證明本研究計算得到的G0是符合實際情況的。

表3 各站點G0與G的線性關系Table 3 The linear relationship between G0 and G of each site

圖3 唐古拉站2015年地表層土壤熱通量（G0）與深度為5 cm的土壤熱通量（G5）散點分布和年均日變化曲線Fig.3 Scatter point distribution and annual mean diurnal variation curve of ground surface soil heat flux（G0）and soil heat flux（G5）with a depth of 5 cm at TGL station in 2015

2.2 遙感反演G0精度驗證

利用計算得到的單點G0值對2015年Ma模型反演得到的估算結果進行精度驗證，采用均方根誤差（root mean square error，RMSE）表征估算的精確程度。計算結果如圖4所示，兩者之間的RMSE為64.81 W·m-2，Ma模型估算值比較接近對應的計算值，決定系數(shù)R2為0.6645，模擬與計算得到的G0之間擬合線的斜率為0.301，出現(xiàn)低估的現(xiàn)象，可能是由于G0與凈輻射之間存在相位差以及輸入遙感數(shù)據(jù)的精度問題所致［21］。

圖4 Ma模型估算地表層土壤熱通量（G0_Ma）與對應地表層土壤熱通量單點計算值（G0_obs）的比較Fig.4 Comparison between ground surface soil heat flux（G0_Ma）estimated by Ma model and corresponding surface soil heat flux single point calculated value（G0_obs）

2.3 地表層土壤熱通量季節(jié)變化分析

2.3.1 單點季節(jié)變化分析 以高原季節(jié)劃分春季（3-5月）、夏季（6-8月）、秋季（9-11月）和冬季（1，2，12月）［39］分析單點G0的季節(jié)變化特征（圖5和表4）。唐古拉與西大灘站（其余站點未列出）2015年表層土壤熱通量G0振幅隨季節(jié)變化（圖5）為夏季較大，冬季較小。這主要是由于G0與凈輻射有極顯著的正相關關系［40］，凈輻射會隨著太陽高度角的變化呈現(xiàn)出顯著的季節(jié)變化特征［8，41］，同時G0會受地表反照率的影響，而地表反照率與植被覆蓋變化呈現(xiàn)相反的變化趨勢［42］。各站點地表層土壤熱通量G0的季節(jié)與年均值如表4所示，G0在春夏兩季均為正值（G0傳播方向以向下為正），說明能量由地上向土壤傳遞，土壤為熱匯；G0在秋冬季除大沙龍站均為負值，說明能量由土壤向地上傳遞，土壤為熱源，此結果與前人研究結果基本一致［43］。從G0年均值可以發(fā)現(xiàn)，各站點有正有負，年均值不超過7 W·m-2。

表4 各站點地表層土壤熱通量的季節(jié)與年均值Table 4 Seasonal and annual average values of ground surface soil heat flux in the surface layer of each station

圖5 唐古拉站、西大灘站地表層土壤熱通量季節(jié)變化Fig.5 Seasonal variation of soil heat flux in surface layer at TGL station and XDT station

2.3.2 G0季節(jié)空間分布特征分析 選擇2015年第113，193，289和353天分別代表春、夏、秋、冬的4 d為高原G0的模擬結果（圖6），G0大小依次為夏季＞春季＞秋季＞冬季。夏季，高原西北地區(qū)的地表層土壤熱通量相對于東南地區(qū)的較高，而冬季則相反。G0的季節(jié)空間分布特征的主要影響因子是植被覆蓋度與地表溫度。夏季，高原的植被覆蓋度分布從西北向東南遞增［44］，且G0因夏季的植被覆蓋度增大而降低［45］。冬季，地表溫度是影響G0的主要因子，高原的地溫由西北向東南逐漸遞增，此空間變化與地表G0有高度一致性。

圖6 2015年第113、193、289和353天地表層土壤熱通量的空間分布Fig.6 Spatial distribution of ground surface soil heat flux on days 113，193，289 and 353 in 2015

2.4 不同草地類型的地表層土壤熱通量

2.4.1 不同草地類型G0年變化 各類草地類型2003-2018年G0的總體均值為40～80 W·m-2，所有草地G0的年變化在2005年最低，均值為55.226 W·m-2，2009年達到最大，均值為63.141 W·m-2，最值之間的差距只有7.915 W·m-2，可見高原草地16年間的G0變化較?。▓D7）。2003-2018年各類草地G0平均值分別為：溫性草原化荒漠類76.557 W·m-2＞溫性草原類72.320 W·m-2＞暖性灌草叢類65.224 W·m-2＞溫性荒漠草原類62.376 W·m-2＞高寒草甸草原類56.034 W·m-2＞沼澤類55.972 W·m-2＞高寒草原類55.623 W·m-2＞高寒荒漠草原類54.226 W·m-2＞山地草甸類50.251 W·m-2＞高寒荒漠類47.207 W·m-2＞高寒草甸類46.118 W·m-2。

圖7 不同草地類型2003-2018年G0均值變化Fig.7 Variation of G0 mean value of different grassland types from 2003 to 2018

總體來看，草地的G0一年內(nèi)呈現(xiàn)出先增后減的變化，峰值基本上集中在5-7月，總體低于140 W·m-2，主要是由于太陽輻射強度與地表溫度的影響（圖8）。一年內(nèi)各類草地類型的平均G0大小依次為溫性草原類71.986 W·m-2，暖性灌草叢類71.425 W·m-2，溫性草原化荒漠類67.833 W·m-2，溫性荒漠草原類63.696 W·m-2，高寒草甸草原類56.202 W·m-2，高寒草原類55.562 W·m-2，沼澤類54.650 W·m-2，高寒荒漠草原類53.121 W·m-2，山地草甸類50.081 W·m-2，高寒荒漠類46.315 W·m-2，高寒草甸類45.444 W·m-2。

圖8 不同草地類型G0的年內(nèi)變化Fig.8 Annual variation of G0 in different grassland types

總的來說，高原各類草地的G0年變化較高的3類草地類型分別是溫性草原化荒漠類、溫性草原類和暖性灌草叢類；G0較低的3類草地是山地草甸類、高寒荒漠類和高寒草甸類。

選擇高原主要的4類草地類型，分別為高寒草甸類、高寒草原類、高寒荒漠草原類和高寒草甸草原類，分析4類草地G0在高原的空間分布情況，如圖9所示。4類草地的G0主要集中在25～100 W·m-2，其中高寒草甸草原類的平均G0最高，高寒草甸類的平均G0最低。高寒草甸主要分布在高原中部偏東，G0沒有明顯的空間分布特征；高寒草原主要分布在高原的西部，G0的空間大小變化從中部向南呈現(xiàn)遞增趨勢；高寒荒漠草原類主要分布在西藏自治區(qū)的北部，G0從東北向西南遞增；高寒草甸草原主要分布在西藏自治區(qū)的中部，G0主要集中于50～75 W·m-2，中間部分較高，四周較低。

圖9 4類草地類型地表層土壤熱通量的空間分布Fig.9 Spatial distribution of ground surface soil heat flux of four grassland types

2.4.2 不同草地類型G0季節(jié)變化 各類草地類型的G0季節(jié)變化呈現(xiàn)出春季到夏季不斷增大，夏季到冬季降低趨勢，且所有草地的G0大小排序為夏季最高，均值93.740 W·m-2，春季次高，均值69.131 W·m-2，秋季較低，均值43.045 W·m-2，冬季最低，均值18.996 W·m-（2圖10）。在春季，溫性草原類G0最高，為88.297 W·m-2，高寒草甸類G0最低，為54.855 W·m-2；在夏季，溫性草原化荒漠類最高，為120.154 W·m-2，高寒草甸類最低，為70.848 W·m-2；在秋季，暖性灌草叢類的G0達到最大，為61.742 W·m-2，高寒荒漠類最低，為28.074 W·m-2；在冬季，暖性灌草叢類的G0最大，為44.155 W·m-2，高寒荒漠類最低，為0.853 W·m-2。

圖10 不同草地類型的G0季節(jié)變化Fig.10 Seasonal variation of G0 in different grassland types

3 討論

本研究通過從不同維度分析高原G0的變化特征，從站點的角度分析，G0均大于不同深度的土壤熱通量值，日變化曲線呈倒“U”形狀，且站點G0季節(jié)振幅變化呈現(xiàn)夏＞春＞秋＞冬，這一結論與前人研究結果一致［38，43］。但是站點年均值G0有正有負且大小在7 W·m-2以下，這與楊成等［43］的研究結果不一致。從空間分布的角度分析，高原G0的季節(jié)變化會受到地表溫度和植被覆蓋度的共同影響，兩因子的影響程度還需要更進一步探討。

從不同時間維度分析高原草地的G0變化。從16年的時間跨度看，每年高原草地G0值變化較小，呈現(xiàn)不規(guī)律的變化趨勢；從一年內(nèi)的變化緯向看，高原草地的整體變化趨勢呈現(xiàn)出先增后降，且峰值出現(xiàn)在5-7月，這與張法偉等［46］發(fā)現(xiàn)高原北部海北高寒矮嵩草（Kobersia humilis）草甸土壤熱通量的年內(nèi)最高值出現(xiàn)在5月有較好的一致性。

全球氣候變化，直接影響大氣-植被-土壤循環(huán)系統(tǒng)，植被生長受到大氣氣候的影響，同時也會受制于土壤水熱的作用，與此同時，植被又會反作用于土壤，所以探討土壤水熱及能量與植被生長的相互關系極具意義。通過探討高原草地G0的變化情況，對研究高原地表能量有所幫助，今后需要繼續(xù)完善對高原草地其他地表熱通量因子（凈輻射通量、感熱通量和潛熱通量）的研究。其次本研究在反演地表層土壤熱通量所選擇的模型與數(shù)據(jù)不是最優(yōu)越的，目前已有國外學者將神經(jīng)網(wǎng)絡模型融入G0的反演［5］，更好地提高G0反演精度。為了更準確地反演高原草地的地表層土壤熱通量，提高高原地表能量平衡模擬精度，為高原草地的可持續(xù)發(fā)展提供更精確可靠的指導，在目前的研究基礎上選擇更精確的模型方法以及加強高原土壤水熱的實測力度將是今后需要努力的方向與目標。

4 結論

本研究基于青藏高原的站點實測數(shù)據(jù)和遙感驅動數(shù)據(jù)進行地表層土壤熱通量的單點計算和Ma模型反演，分別從點和面尺度分析了高原地表層土壤熱通量時空分布特征變化，并且分析了不同草地類型的G0變化，得出以下主要結論：

1）站點地表層土壤熱通量比不同深度的土壤熱通量值大。G0的日變化曲線呈倒“U”形狀，在夜晚相較于白天變化較為平緩。

2）站點地表層土壤熱通量的季節(jié)振幅變化呈現(xiàn)夏＞春＞秋＞冬。春夏季G0均值為正值，土壤為熱匯；秋冬季G0均值基本為負值，土壤為熱源。從空間分布分析，夏季，高原西北地區(qū)的地表層土壤熱通量相對于東南地區(qū)的較高，而冬季則相反。

3）高原草地的土壤熱通量值為40～80 W·m-2，16年各類草地G0平均值最高的是溫性草原化荒漠類（76.557 W·m-2），最低的是高寒草甸類（46.118 W·m-2）。

4）高原草地的G0一年內(nèi)呈現(xiàn)出先增后減的變化趨勢，總體低于140 W·m-2。高原草地的G0有明顯的季節(jié)變化?？傮w季節(jié)變化規(guī)律呈現(xiàn)夏＞春＞秋＞冬，且均值大小分別為93.740，69.131，43.045，18.996 W·m-2。