李穎,吳靜,李純斌,秦格霞
(甘肅農(nóng)業(yè)大學(xué)資源與環(huán)境學(xué)院,甘肅 蘭州 730070)
青藏高原(以下簡(jiǎn)稱高原)占中國(guó)陸地面積的1/4,海拔均在4000 m以上,是世界上海拔最高面積最大的高原[1]。在全球氣候變暖的大背景下,高原的地表能量研究成為熱點(diǎn)課題[2-4]。地表層土壤熱通量(ground surface soil heat flux,G0)是地表能量平衡(surface energy balance,SEB)的重要參數(shù)之一,它可以解釋地表淺層與深層之間傳輸?shù)哪芰浚?]。同時(shí)對(duì)地氣間的物質(zhì)能量傳輸及分配有重要影響[6]。此外,G0也是地表蒸散發(fā)模擬的重要因子,會(huì)極大地影響地表蒸散發(fā)的模擬精度[7]。一直以來(lái)對(duì)高原G0的研究從未中斷,高原腹地的G0主要受到地氣溫差的影響,同時(shí)下墊面狀況也起到重要作用[8];利用超高分辨率掃描輻射計(jì)(advanced very high resolution radiometer,AVHRR)和中分辨率成像光譜儀(moderate resolution image spectroradiometer,MODIS)反演地表特征參數(shù)和地表通量,發(fā)現(xiàn)G0空間分布與地表特征參數(shù)(植被指數(shù)、植被覆蓋度和地表反照率)的空間分布有很好的對(duì)應(yīng)關(guān)系[9]。由此可見(jiàn),高原G0與下墊面有非常密切的關(guān)系。高原分布著占全國(guó)草地總面積38%的天然草地,不僅具有防風(fēng)固沙、涵養(yǎng)水源、調(diào)節(jié)碳循環(huán)及氣候變化、維護(hù)生物多樣性等多種生態(tài)服務(wù)功能[10],同時(shí)還是當(dāng)?shù)鼐用裆姘l(fā)展和畜牧業(yè)發(fā)展的基礎(chǔ)[11]。通過(guò)研究高原草地G0,不僅有助于加深對(duì)高原G0的理解,而且對(duì)高原土壤-植被-大氣系統(tǒng)的能量平衡研究有所幫助。
目前對(duì)高原地表層土壤熱通量G0的研究主要有兩種手段,一是基于站點(diǎn)觀測(cè)數(shù)據(jù)估算分析G0,二是結(jié)合遙感影像數(shù)據(jù)分析時(shí)空尺度較大的G0,中心思想是通過(guò)獲取G0與凈輻射(net radiation,Rn)的比率,以此估計(jì)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]。有學(xué)者對(duì)這些模型做了對(duì)比分析,對(duì)于特殊的高原氣候條件,反演地表層土壤熱通量最適用的方法是Ma模型[15]。由于高原的地理氣候條件,實(shí)測(cè)數(shù)據(jù)難以獲取,因此利用遙感數(shù)據(jù)估算G0是非常有效的手段。本研究利用Ma模型反演高原2003-2018年地表層土壤熱通量(G0),分析2003-2018年高原G0的時(shí)空變化特征,探討不同草地類型的G0年際變化以及季節(jié)變化。
青藏高原(26°00′12″-39°46′50″N,73°18′52″-104°46′59″E)總體地勢(shì)西高東低,地形復(fù)雜,地勢(shì)險(xiǎn)峻多變。草地是高原主要的生態(tài)系統(tǒng)類型,分布情況如圖1所示,主要以高寒草原類和高寒草甸類為主,占高原草地面積的73.8%。
1.2.1 草地類型數(shù)據(jù) 草地類型數(shù)據(jù)來(lái)自蘇大學(xué)[16]團(tuán)隊(duì)編制的1∶1000000的中國(guó)草地類型數(shù)據(jù),利用高原行政區(qū)域提取高原的草地類型,分類合并,得到高原的草地類型圖(圖1)。
圖1 青藏高原草地類型及觀測(cè)點(diǎn)分布Fig.1 Grassland types and distribution of observation points on the Qinghai Tibet Plateau
1.2.2 站點(diǎn)數(shù)據(jù) 選擇7個(gè)高原觀測(cè)土壤熱通量的站點(diǎn)(圖1和表1),其中大沙龍站,阿柔站,峨堡站,景陽(yáng)嶺站,埡口站數(shù)據(jù)來(lái)源于國(guó)家青藏高原科學(xué)數(shù)據(jù)中心(http://data.tpdc.ac.cn/zh-hans/),本數(shù)據(jù)由“黑河生態(tài)水文遙感試驗(yàn)(HiWATER)”產(chǎn)生[17-18];西大灘站,唐古拉站數(shù)據(jù)來(lái)源于國(guó)家冰川凍土沙漠科學(xué)數(shù)據(jù)中心(http://www.crensed.ac.cn/portal/)[19-20]。唐古拉站和西大灘站的地表溫度(temperature,T)用史蒂芬-玻爾茲曼公式計(jì)算獲得[21]。
表1 觀測(cè)站點(diǎn)信息Table 1 The information of observation sites
式中:T為地表溫度(單位:K);Rul和Rdl分別為地表向上長(zhǎng)波輻射和地表向下長(zhǎng)波輻射(單位:W·m-2);εg為地表比輻射率,取值為0.96;σ為Stefan-Boltzmann常數(shù),取5.67×10-8W·m-2·K-4。
1.2.3 遙感數(shù)據(jù) 地表層土壤熱通量的遙感反演數(shù)據(jù)包括中國(guó)西部1 km全天候地表溫度數(shù)據(jù)V1,覆蓋整個(gè)青藏高原,時(shí)間分辨率為逐日,空間分辨率為1 km[22];MODIS數(shù)據(jù)來(lái)自于美國(guó)航空航天局,包括MOD09GA和MOD13Q1。其中MOD09GA提供1~7波段的表面光譜反射率估計(jì)值,主要用于地表反照率的計(jì)算;MOD13Q1提供NDVI、紅光和近紅外波段地表反射率值,NDVI用于地表發(fā)射率的反演,紅光和近紅外波段的地表反射率用于反演區(qū)域改良的土壤調(diào)整植被指數(shù)(modified soil adjusted vegetation index,MSAVI);中國(guó)區(qū)域地面氣象要素驅(qū)動(dòng)數(shù)據(jù)集(1979-2018)來(lái)自國(guó)家青藏高原科學(xué)數(shù)據(jù)中心,該數(shù)據(jù)集主要提供近地面氣溫、近地面氣壓、近地面空氣比濕、近地面全風(fēng)速、地面向下短波輻射、地面向下長(zhǎng)波輻射和地面降水率共7個(gè)要素[23-26]。所需遙感數(shù)據(jù)具體介紹如表2所示。
表2 遙感數(shù)據(jù)信息Table 2 The information of remote sensing data
首先利用單點(diǎn)實(shí)測(cè)數(shù)據(jù),計(jì)算單點(diǎn)地表土壤熱通量;其次根據(jù)遙感反演模型,利用遙感影像數(shù)據(jù)反演得到高原地表層土壤熱通量,分析不同草地類型的地表層土壤熱通量,具體流程如圖2所示。
圖2 技術(shù)路線Fig.2 Technology road
1.3.1 單點(diǎn)實(shí)測(cè)數(shù)據(jù)計(jì)算方法 采用以下公式計(jì)算站點(diǎn)地表層土壤熱通量:
式中:t為時(shí)間(單位:s);z為土壤深度(向下為正,單位:m);T為土壤溫度(單位:K);Cs為土壤比熱容(單位:J·kg-1·K-1);ρs為土壤密度(單位:kg·m-3);λs為土壤熱傳導(dǎo)系數(shù)(單位:W·K-1·m-1);G為土壤熱通量(向下為正,單位:W·m-2)[27]。
結(jié)合公式2)與3),且考慮青藏高原近地表存在土壤凍融現(xiàn)象,因此利用5 cm的土壤熱通量計(jì)算地表層土壤熱通量(G0)[28]:
式中:zref為參考深度(基于本研究的實(shí)測(cè)數(shù)據(jù),不同站點(diǎn)選擇不同的參考深度);G(zref)為參考深度處的土壤熱通量;ρw為液態(tài)水的密度,為1.0×103kg·m-3;Lf為水的凍結(jié)-融化潛熱常數(shù),為3.337×105J·kg-1;θw為體積含水量;Tav為5 cm土壤溫度與地表溫度的平均值。在土壤的整個(gè)凍融過(guò)程中,對(duì)不同凍融階段進(jìn)行不同的計(jì)算[29]:
1)當(dāng)5 cm土壤層處于完全融化階段時(shí),土壤體積熱容量計(jì)算公式如下:
式中:ρ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)當(dāng)5 cm深度以上的土壤層處于凍融循環(huán)階段時(shí),此時(shí)土壤體積熱容量為:
式中:ρiCi為冰的熱容量,約為1.89×106J·m-3·K-1;θi5cm為5 cm處的土壤含冰量。在凍融循環(huán)階段,土壤中的含冰量可以根據(jù)對(duì)應(yīng)深度的土壤中相鄰時(shí)刻之間未凍水含量來(lái)計(jì)算[31],公式如下:
3)當(dāng)土層5 cm以上完全凍結(jié)時(shí),可以將土壤含水量和含冰量視為常數(shù)[32]。土壤中的含冰量由凍結(jié)開(kāi)始前和結(jié)束后含水量的差值求得,其土壤體積熱容量仍用公式6)進(jìn)行計(jì)算。
1.3.2 遙感反演方法 Ma總模型公式如下:
式中:Ts為地表溫度(單位:℃);α是地表反照率;MSAVI為改良的土壤調(diào)整植被指數(shù)[15];Rn是凈輻射(單位:W·m-2)。
地表反照率α基于MOD09GA的波段反射率數(shù)據(jù)計(jì)算得出[33]:
式中:R1、R2、R3、R4、R5和R7為MOD09GA的第1、2、3、4、5和7波段的地表反射率。
MSAVI通過(guò)以下公式計(jì)算得出[34]:
式中:RED和NIR分別為MOD13Q1的紅光波段和近紅外波段的地表反射率值。
凈輻射Rn可以通過(guò)以下公式計(jì)算獲取[35]:
式中:向下短波輻射(downward short wave radiation,DSR)和向下長(zhǎng)波輻射(downward long wave radiation,DLR)分別是向下短波輻射和向下長(zhǎng)波輻射(單位:W·m-2);εs是地表發(fā)射率,根據(jù)以下公式計(jì)算[36]:
式中:NDVI為歸一化植被指數(shù),通過(guò)植被指數(shù)產(chǎn)品MOD13Q1提?。籪c為植被蓋度,可以根據(jù)NDVI估算,公式如下[37]:
式中:NDVImax和NDVImin分別是完全植被覆蓋和裸露土壤的NDVI值。
由于G0無(wú)法直接觀測(cè)得到,因此本研究基于站點(diǎn)淺層的土壤熱通量觀測(cè)值G檢驗(yàn)站點(diǎn)G0的計(jì)算結(jié)果。圖3a是唐古拉站點(diǎn)2015年G0與5 cm處的土壤熱通量G5之間的散點(diǎn)分布關(guān)系,各站點(diǎn)G0與淺層土壤熱通量G的線性關(guān)系式以及R2的值如表3所示,除阿柔站,其余站點(diǎn)的G0與G無(wú)截距線性關(guān)系,斜率略大于1,說(shuō)明G0略大于G5,兩者之間的R2均大于0.8。圖3b是唐古拉站點(diǎn)G0與G52015年年均日變化曲線,兩者整體變化曲線呈現(xiàn)倒“U”形狀,在夜晚相較于白天變化較為平緩。在7:00之前,G5略大于G0,并且G0比G5首先達(dá)到最大值。以上結(jié)論與之前的學(xué)者研究結(jié)果一致[38],可以證明本研究計(jì)算得到的G0是符合實(shí)際情況的。
表3 各站點(diǎn)G0與G的線性關(guān)系Table 3 The linear relationship between G0 and G of each site
圖3 唐古拉站2015年地表層土壤熱通量(G0)與深度為5 cm的土壤熱通量(G5)散點(diǎn)分布和年均日變化曲線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
利用計(jì)算得到的單點(diǎn)G0值對(duì)2015年Ma模型反演得到的估算結(jié)果進(jìn)行精度驗(yàn)證,采用均方根誤差(root mean square error,RMSE)表征估算的精確程度。計(jì)算結(jié)果如圖4所示,兩者之間的RMSE為64.81 W·m-2,Ma模型估算值比較接近對(duì)應(yīng)的計(jì)算值,決定系數(shù)R2為0.6645,模擬與計(jì)算得到的G0之間擬合線的斜率為0.301,出現(xiàn)低估的現(xiàn)象,可能是由于G0與凈輻射之間存在相位差以及輸入遙感數(shù)據(jù)的精度問(wèn)題所致[21]。
圖4 Ma模型估算地表層土壤熱通量(G0_Ma)與對(duì)應(yīng)地表層土壤熱通量單點(diǎn)計(jì)算值(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.1 單點(diǎn)季節(jié)變化分析 以高原季節(jié)劃分春季(3-5月)、夏季(6-8月)、秋季(9-11月)和冬季(1,2,12月)[39]分析單點(diǎn)G0的季節(jié)變化特征(圖5和表4)。唐古拉與西大灘站(其余站點(diǎn)未列出)2015年表層土壤熱通量G0振幅隨季節(jié)變化(圖5)為夏季較大,冬季較小。這主要是由于G0與凈輻射有極顯著的正相關(guān)關(guān)系[40],凈輻射會(huì)隨著太陽(yáng)高度角的變化呈現(xiàn)出顯著的季節(jié)變化特征[8,41],同時(shí)G0會(huì)受地表反照率的影響,而地表反照率與植被覆蓋變化呈現(xiàn)相反的變化趨勢(shì)[42]。各站點(diǎn)地表層土壤熱通量G0的季節(jié)與年均值如表4所示,G0在春夏兩季均為正值(G0傳播方向以向下為正),說(shuō)明能量由地上向土壤傳遞,土壤為熱匯;G0在秋冬季除大沙龍站均為負(fù)值,說(shuō)明能量由土壤向地上傳遞,土壤為熱源,此結(jié)果與前人研究結(jié)果基本一致[43]。從G0年均值可以發(fā)現(xiàn),各站點(diǎn)有正有負(fù),年均值不超過(guò)7 W·m-2。
表4 各站點(diǎn)地表層土壤熱通量的季節(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的模擬結(jié)果(圖6),G0大小依次為夏季>春季>秋季>冬季。夏季,高原西北地區(qū)的地表層土壤熱通量相對(duì)于東南地區(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.1 不同草地類型G0年變化 各類草地類型2003-2018年G0的總體均值為40~80 W·m-2,所有草地G0的年變化在2005年最低,均值為55.226 W·m-2,2009年達(dá)到最大,均值為63.141 W·m-2,最值之間的差距只有7.915 W·m-2,可見(jiàn)高原草地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
總體來(lái)看,草地的G0一年內(nèi)呈現(xiàn)出先增后減的變化,峰值基本上集中在5-7月,總體低于140 W·m-2,主要是由于太陽(yáng)輻射強(qiáng)度與地表溫度的影響(圖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
總的來(lái)說(shuō),高原各類草地的G0年變化較高的3類草地類型分別是溫性草原化荒漠類、溫性草原類和暖性灌草叢類;G0較低的3類草地是山地草甸類、高寒荒漠類和高寒草甸類。
選擇高原主要的4類草地類型,分別為高寒草甸類、高寒草原類、高寒荒漠草原類和高寒草甸草原類,分析4類草地G0在高原的空間分布情況,如圖9所示。4類草地的G0主要集中在25~100 W·m-2,其中高寒草甸草原類的平均G0最高,高寒草甸類的平均G0最低。高寒草甸主要分布在高原中部偏東,G0沒(méi)有明顯的空間分布特征;高寒草原主要分布在高原的西部,G0的空間大小變化從中部向南呈現(xiàn)遞增趨勢(shì);高寒荒漠草原類主要分布在西藏自治區(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)出春季到夏季不斷增大,夏季到冬季降低趨勢(shì),且所有草地的G0大小排序?yàn)橄募咀罡撸?3.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達(dá)到最大,為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
本研究通過(guò)從不同維度分析高原G0的變化特征,從站點(diǎn)的角度分析,G0均大于不同深度的土壤熱通量值,日變化曲線呈倒“U”形狀,且站點(diǎn)G0季節(jié)振幅變化呈現(xiàn)夏>春>秋>冬,這一結(jié)論與前人研究結(jié)果一致[38,43]。但是站點(diǎn)年均值G0有正有負(fù)且大小在7 W·m-2以下,這與楊成等[43]的研究結(jié)果不一致。從空間分布的角度分析,高原G0的季節(jié)變化會(huì)受到地表溫度和植被覆蓋度的共同影響,兩因子的影響程度還需要更進(jìn)一步探討。
從不同時(shí)間維度分析高原草地的G0變化。從16年的時(shí)間跨度看,每年高原草地G0值變化較小,呈現(xiàn)不規(guī)律的變化趨勢(shì);從一年內(nèi)的變化緯向看,高原草地的整體變化趨勢(shì)呈現(xiàn)出先增后降,且峰值出現(xiàn)在5-7月,這與張法偉等[46]發(fā)現(xiàn)高原北部海北高寒矮嵩草(Kobersia humilis)草甸土壤熱通量的年內(nèi)最高值出現(xiàn)在5月有較好的一致性。
全球氣候變化,直接影響大氣-植被-土壤循環(huán)系統(tǒng),植被生長(zhǎng)受到大氣氣候的影響,同時(shí)也會(huì)受制于土壤水熱的作用,與此同時(shí),植被又會(huì)反作用于土壤,所以探討土壤水熱及能量與植被生長(zhǎng)的相互關(guān)系極具意義。通過(guò)探討高原草地G0的變化情況,對(duì)研究高原地表能量有所幫助,今后需要繼續(xù)完善對(duì)高原草地其他地表熱通量因子(凈輻射通量、感熱通量和潛熱通量)的研究。其次本研究在反演地表層土壤熱通量所選擇的模型與數(shù)據(jù)不是最優(yōu)越的,目前已有國(guó)外學(xué)者將神經(jīng)網(wǎng)絡(luò)模型融入G0的反演[5],更好地提高G0反演精度。為了更準(zhǔn)確地反演高原草地的地表層土壤熱通量,提高高原地表能量平衡模擬精度,為高原草地的可持續(xù)發(fā)展提供更精確可靠的指導(dǎo),在目前的研究基礎(chǔ)上選擇更精確的模型方法以及加強(qiáng)高原土壤水熱的實(shí)測(cè)力度將是今后需要努力的方向與目標(biāo)。
本研究基于青藏高原的站點(diǎn)實(shí)測(cè)數(shù)據(jù)和遙感驅(qū)動(dòng)數(shù)據(jù)進(jìn)行地表層土壤熱通量的單點(diǎn)計(jì)算和Ma模型反演,分別從點(diǎn)和面尺度分析了高原地表層土壤熱通量時(shí)空分布特征變化,并且分析了不同草地類型的G0變化,得出以下主要結(jié)論:
1)站點(diǎn)地表層土壤熱通量比不同深度的土壤熱通量值大。G0的日變化曲線呈倒“U”形狀,在夜晚相較于白天變化較為平緩。
2)站點(diǎn)地表層土壤熱通量的季節(jié)振幅變化呈現(xiàn)夏>春>秋>冬。春夏季G0均值為正值,土壤為熱匯;秋冬季G0均值基本為負(fù)值,土壤為熱源。從空間分布分析,夏季,高原西北地區(qū)的地表層土壤熱通量相對(duì)于東南地區(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)出先增后減的變化趨勢(shì),總體低于140 W·m-2。高原草地的G0有明顯的季節(jié)變化??傮w季節(jié)變化規(guī)律呈現(xiàn)夏>春>秋>冬,且均值大小分別為93.740,69.131,43.045,18.996 W·m-2。