2003-2018年青藏高原草地的地表层土壤热通量时空变化

李颖，吴静，李纯斌，秦格霞

（甘肃农业大学资源与环境学院，甘肃 兰州 730070）

青藏高原（以下简称高原）占中国陆地面积的1/4，海拔均在4000 m以上，是世界上海拔最高面积最大的高原［1］。在全球气候变暖的大背景下，高原的地表能量研究成为热点课题［2-4］。地表层土壤热通量（ground surface soil heat flux，G0）是地表能量平衡（surface energy balance，SEB）的重要参数之一，它可以解释地表浅层与深层之间传输的能量［5］。同时对地气间的物质能量传输及分配有重要影响［6］。此外，G0也是地表蒸散发模拟的重要因子，会极大地影响地表蒸散发的模拟精度［7］。一直以来对高原G0的研究从未中断，高原腹地的G0主要受到地气温差的影响，同时下垫面状况也起到重要作用［8］；利用超高分辨率扫描辐射计（advanced very high resolution radiometer，AVHRR）和中分辨率成像光谱仪（moderate resolution image spectroradiometer，MODIS）反演地表特征参数和地表通量，发现G0空间分布与地表特征参数（植被指数、植被覆盖度和地表反照率）的空间分布有很好的对应关系［9］。由此可见，高原G0与下垫面有非常密切的关系。高原分布着占全国草地总面积38%的天然草地，不仅具有防风固沙、涵养水源、调节碳循环及气候变化、维护生物多样性等多种生态服务功能［10］，同时还是当地居民生存发展和畜牧业发展的基础［11］。通过研究高原草地G0，不仅有助于加深对高原G0的理解，而且对高原土壤-植被-大气系统的能量平衡研究有所帮助。

目前对高原地表层土壤热通量G0的研究主要有两种手段，一是基于站点观测数据估算分析G0，二是结合遥感影像数据分析时空尺度较大的G0，中心思想是通过获取G0与净辐射（net radiation，Rn）的比率，以此估计G0。利用不同的地表参数估算G0/Rn，如基于植被参数的模型，有归一化植被指数（normalized difference vegetation index，NDVI）和叶面积指数（leaf area index，LAI）模型［12］；除此之外考虑地表温度和地表反照率的模型，如地表能量平衡模型（surface energy balance model，SEBAL）［13］和Ma模型［14］。有学者对这些模型做了对比分析，对于特殊的高原气候条件，反演地表层土壤热通量最适用的方法是Ma模型［15］。由于高原的地理气候条件，实测数据难以获取，因此利用遥感数据估算G0是非常有效的手段。本研究利用Ma模型反演高原2003-2018年地表层土壤热通量（G0），分析2003-2018年高原G0的时空变化特征，探讨不同草地类型的G0年际变化以及季节变化。

1 材料与方法

1.1 研究区概况

青藏高原（26°00′12″-39°46′50″N，73°18′52″-104°46′59″E）总体地势西高东低，地形复杂，地势险峻多变。草地是高原主要的生态系统类型，分布情况如图1所示，主要以高寒草原类和高寒草甸类为主，占高原草地面积的73.8%。

1.2 数据来源及处理

1.2.1 草地类型数据 草地类型数据来自苏大学［16］团队编制的1∶1000000的中国草地类型数据，利用高原行政区域提取高原的草地类型，分类合并，得到高原的草地类型图（图1）。

图1 青藏高原草地类型及观测点分布Fig.1 Grassland types and distribution of observation points on the Qinghai Tibet Plateau

1.2.2 站点数据 选择7个高原观测土壤热通量的站点（图1和表1），其中大沙龙站，阿柔站，峨堡站，景阳岭站，垭口站数据来源于国家青藏高原科学数据中心（http：//data.tpdc.ac.cn/zh-hans/），本数据由“黑河生态水文遥感试验（HiWATER）”产生［17-18］；西大滩站，唐古拉站数据来源于国家冰川冻土沙漠科学数据中心（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常数，取5.67×10-8W·m-2·K-4。

1.2.3 遥感数据 地表层土壤热通量的遥感反演数据包括中国西部1 km全天候地表温度数据V1，覆盖整个青藏高原，时间分辨率为逐日，空间分辨率为1 km［22］；MODIS数据来自于美国航空航天局，包括MOD09GA和MOD13Q1。其中MOD09GA提供1～7波段的表面光谱反射率估计值，主要用于地表反照率的计算；MOD13Q1提供NDVI、红光和近红外波段地表反射率值，NDVI用于地表发射率的反演，红光和近红外波段的地表反射率用于反演区域改良的土壤调整植被指数（modified soil adjusted vegetation index，MSAVI）；中国区域地面气象要素驱动数据集（1979-2018）来自国家青藏高原科学数据中心，该数据集主要提供近地面气温、近地面气压、近地面空气比湿、近地面全风速、地面向下短波辐射、地面向下长波辐射和地面降水率共7个要素［23-26］。所需遥感数据具体介绍如表2所示。

表2 遥感数据信息Table 2 The information of remote sensing data

1.3 研究方法

首先利用单点实测数据，计算单点地表土壤热通量；其次根据遥感反演模型，利用遥感影像数据反演得到高原地表层土壤热通量，分析不同草地类型的地表层土壤热通量，具体流程如图2所示。

图2 技术路线Fig.2 Technology road

1.3.1 单点实测数据计算方法 采用以下公式计算站点地表层土壤热通量：

式中：t为时间（单位：s）；z为土壤深度（向下为正，单位：m）；T为土壤温度（单位：K）；Cs为土壤比热容（单位：J·kg-1·K-1）；ρs为土壤密度（单位：kg·m-3）；λs为土壤热传导系数（单位：W·K-1·m-1）；G为土壤热通量（向下为正，单位：W·m-2）［27］。

结合公式2）与3），且考虑青藏高原近地表存在土壤冻融现象，因此利用5 cm的土壤热通量计算地表层土壤热通量（G0）［28］：

式中：zref为参考深度（基于本研究的实测数据，不同站点选择不同的参考深度）；G(zref)为参考深度处的土壤热通量；ρw为液态水的密度，为1.0×103kg·m-3；Lf为水的冻结-融化潜热常数，为3.337×105J·kg-1；θw为体积含水量；Tav为5 cm土壤温度与地表温度的平均值。在土壤的整个冻融过程中，对不同冻融阶段进行不同的计算［29］：

1）当5 cm土壤层处于完全融化阶段时，土壤体积热容量计算公式如下：

式中：ρdryCdry为干土的热容量，约为0.90×106J·m-3·K-1［30］；ρwCw为液态水的热容量，约为4.2×106J·m-3·K-1；θw5cm为深度5 cm处未冻水体积含水量（单位：m3·m-3）。

2）当5 cm深度以上的土壤层处于冻融循环阶段时，此时土壤体积热容量为：

式中：ρiCi为冰的热容量，约为1.89×106J·m-3·K-1；θi5cm为5 cm处的土壤含冰量。在冻融循环阶段，土壤中的含冰量可以根据对应深度的土壤中相邻时刻之间未冻水含量来计算［31］，公式如下：

3）当土层5 cm以上完全冻结时，可以将土壤含水量和含冰量视为常数［32］。土壤中的含冰量由冻结开始前和结束后含水量的差值求得，其土壤体积热容量仍用公式6）进行计算。

1.3.2 遥感反演方法 Ma总模型公式如下：

式中：Ts为地表温度（单位：℃）；α是地表反照率；MSAVI为改良的土壤调整植被指数［15］；Rn是净辐射（单位：W·m-2）。

地表反照率α基于MOD09GA的波段反射率数据计算得出［33］：

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

MSAVI通过以下公式计算得出［34］：

式中：RED和NIR分别为MOD13Q1的红光波段和近红外波段的地表反射率值。

净辐射Rn可以通过以下公式计算获取［35］：

式中：向下短波辐射（downward short wave radiation，DSR）和向下长波辐射（downward long wave radiation，DLR）分别是向下短波辐射和向下长波辐射（单位：W·m-2）；εs是地表发射率，根据以下公式计算［36］：

式中：NDVI为归一化植被指数，通过植被指数产品MOD13Q1提取；fc为植被盖度，可以根据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年年均日变化曲线，两者整体变化曲线呈现倒“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模型估算值比较接近对应的计算值，决定系数R2为0.6645，模拟与计算得到的G0之间拟合线的斜率为0.301，出现低估的现象，可能是由于G0与净辐射之间存在相位差以及输入遥感数据的精度问题所致［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 地表层土壤热通量季节变化分析

2.3.1 单点季节变化分析 以高原季节划分春季（3-5月）、夏季（6-8月）、秋季（9-11月）和冬季（1，2，12月）［39］分析单点G0的季节变化特征（图5和表4）。唐古拉与西大滩站（其余站点未列出）2015年表层土壤热通量G0振幅随季节变化（图5）为夏季较大，冬季较小。这主要是由于G0与净辐射有极显著的正相关关系［40］，净辐射会随着太阳高度角的变化呈现出显著的季节变化特征［8，41］，同时G0会受地表反照率的影响，而地表反照率与植被覆盖变化呈现相反的变化趋势［42］。各站点地表层土壤热通量G0的季节与年均值如表4所示，G0在春夏两季均为正值（G0传播方向以向下为正），说明能量由地上向土壤传递，土壤为热汇；G0在秋冬季除大沙龙站均为负值，说明能量由土壤向地上传递，土壤为热源，此结果与前人研究结果基本一致［43］。从G0年均值可以发现，各站点有正有负，年均值不超过7 W·m-2。

表4 各站点地表层土壤热通量的季节与年均值Table 4 Seasonal and annual average values of ground surface soil heat flux in the surface layer of each station

图5 唐古拉站、西大滩站地表层土壤热通量季节变化Fig.5 Seasonal variation of soil heat flux in surface layer at TGL station and XDT station

2.3.2 G0季节空间分布特征分析 选择2015年第113，193，289和353天分别代表春、夏、秋、冬的4 d为高原G0的模拟结果（图6），G0大小依次为夏季＞春季＞秋季＞冬季。夏季，高原西北地区的地表层土壤热通量相对于东南地区的较高，而冬季则相反。G0的季节空间分布特征的主要影响因子是植被覆盖度与地表温度。夏季，高原的植被覆盖度分布从西北向东南递增［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变化较小（图7）。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一年内呈现出先增后减的变化，峰值基本上集中在5-7月，总体低于140 W·m-2，主要是由于太阳辐射强度与地表温度的影响（图8）。一年内各类草地类型的平均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的年内变化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的空间大小变化从中部向南呈现递增趋势；高寒荒漠草原类主要分布在西藏自治区的北部，G0从东北向西南递增；高寒草甸草原主要分布在西藏自治区的中部，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季节变化 各类草地类型的G0季节变化呈现出春季到夏季不断增大，夏季到冬季降低趋势，且所有草地的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季节变化Fig.10 Seasonal variation of G0 in different grassland types

3 讨论

本研究通过从不同维度分析高原G0的变化特征，从站点的角度分析，G0均大于不同深度的土壤热通量值，日变化曲线呈倒“U”形状，且站点G0季节振幅变化呈现夏＞春＞秋＞冬，这一结论与前人研究结果一致［38，43］。但是站点年均值G0有正有负且大小在7 W·m-2以下，这与杨成等［43］的研究结果不一致。从空间分布的角度分析，高原G0的季节变化会受到地表温度和植被覆盖度的共同影响，两因子的影响程度还需要更进一步探讨。

从不同时间维度分析高原草地的G0变化。从16年的时间跨度看，每年高原草地G0值变化较小，呈现不规律的变化趋势；从一年内的变化纬向看，高原草地的整体变化趋势呈现出先增后降，且峰值出现在5-7月，这与张法伟等［46］发现高原北部海北高寒矮嵩草（Kobersia humilis）草甸土壤热通量的年内最高值出现在5月有较好的一致性。

全球气候变化，直接影响大气-植被-土壤循环系统，植被生长受到大气气候的影响，同时也会受制于土壤水热的作用，与此同时，植被又会反作用于土壤，所以探讨土壤水热及能量与植被生长的相互关系极具意义。通过探讨高原草地G0的变化情况，对研究高原地表能量有所帮助，今后需要继续完善对高原草地其他地表热通量因子（净辐射通量、感热通量和潜热通量）的研究。其次本研究在反演地表层土壤热通量所选择的模型与数据不是最优越的，目前已有国外学者将神经网络模型融入G0的反演［5］，更好地提高G0反演精度。为了更准确地反演高原草地的地表层土壤热通量，提高高原地表能量平衡模拟精度，为高原草地的可持续发展提供更精确可靠的指导，在目前的研究基础上选择更精确的模型方法以及加强高原土壤水热的实测力度将是今后需要努力的方向与目标。

4 结论

本研究基于青藏高原的站点实测数据和遥感驱动数据进行地表层土壤热通量的单点计算和Ma模型反演，分别从点和面尺度分析了高原地表层土壤热通量时空分布特征变化，并且分析了不同草地类型的G0变化，得出以下主要结论：

1）站点地表层土壤热通量比不同深度的土壤热通量值大。G0的日变化曲线呈倒“U”形状，在夜晚相较于白天变化较为平缓。

2）站点地表层土壤热通量的季节振幅变化呈现夏＞春＞秋＞冬。春夏季G0均值为正值，土壤为热汇；秋冬季G0均值基本为负值，土壤为热源。从空间分布分析，夏季，高原西北地区的地表层土壤热通量相对于东南地区的较高，而冬季则相反。

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

4）高原草地的G0一年内呈现出先增后减的变化趋势，总体低于140 W·m-2。高原草地的G0有明显的季节变化。总体季节变化规律呈现夏＞春＞秋＞冬，且均值大小分别为93.740，69.131，43.045，18.996 W·m-2。